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ON FINITE ELASTIC DEFORMATIONS WITH A 
PERTURBED STRAIN-ENERGY FUNCTIONt 
By A. J. M. SPENCERt (Brown University) 
[Received 20 February 1958] 


SUMMARY 

In order to obtain analytical solutions to problems involving finite deformations 
of elastic solids it is often necessary to adopt a specific form for the strain-energy 
function. In this paper it is supposed that such a strain-energy function W is 
perturbed to a new strain-energy function W + «W’, with the result that an additional 
small deformation is superimposed on the existing finite deformation. Stress-strain 
relations and equations of equilibrium are formulated which, with suitable boundary 
conditions and a knowledge of the original finite deformation, make possible the 
determination of the displacements and stress-components that are superimposed 
when W is perturbed. The theory is illustrated by application to the problem of 
simultaneous extension, inflation, and shear of a cylindrical tube, and also to the 
special case in which the tube undergoes shear only. In these applications W is 
supposed to have the Mooney form, while W’ is an arbitrary function of the strain 
invariants. 


1. Introduction 

In recent years a number of solutions have been obtained to problems 
involving finite deformations of elastic solids. In some particular problems 
it is possible to obtain complete solutions determining the deformation 
and state of stress in a body without specifying any particular form for 
the strain-energy function of the elastic solid. In other cases, the equations 
of equilibrium are such that they can only be integrated if a specific form 
of the strain-energy function is assumed. This situation arises, for 
example, in the simultaneous extension, inflation, and shear of a cylindrical 
tube. In such cases, the forms most often adopted for the strain-energy 
function W are the form postulated by Mooney (1), 


W = C,(1,—3)+C,(1,—3), (1.1) 
and the neo-Hookean form 
W = C,(I,—3). (1.2) 
Here J, and J, are invariants of the strain tensor to be defined later. 
C, and ©, are positive constants. 
Experiments on vulcanized rubbers by Rivlin and Saunders (2) and by 
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Gent and Rivlin (3) show that for these materials the strain-energy func- 
tion departs slightly from the form (1.1) and for moderately large deforma- 
tions may be expressed by a formula of the type 


W = C(,—3)+f(h—3), (1.3) 
where f is a decreasing function of (/,—3). It is therefore of interest to 
investigate the theory of the deformation of elastic materials with a strain- 
energy function of the form 

W* = W+eW’, (1.4) 
where W is a strain-energy function for which, in a given problem, an 
explicit solution can be found, and «W’ represents a small perturbing 
strain-energy function. By studying this form we may to some extent 
remove the restrictions imposed when the exact solution of the equations 


requires the assumption of some special form for the strain-energy func- 


tion. With a suitable choice of W and W’, the theory we shall develop 
may also be applicable to the behaviour of certain real materials. In 
applications, W will usually have the form (1.1) or (1.2), but the theory 
does not restrict W to any particular form. 

A brief summary of the notation and equations of the theory of finite 
elastic deformations is presented in section 2. This is followed, in section 
3, by a consideration of the effect of replacing the strain-energy function 
W of a deformed elastic body B by a new strain-energy function W+«W’. 
This must, in general, result in a further deformation of the body into 
a second deformed body B’. It is supposed that the deformation and the 
state of stress of the body B is completely determined, and we assume 
that the additional displacements that occur when B is deformed to B’ 
are everywhere of order «. The theory then closely resembles the theory 
of small elastic deformations superposed on finite elastic deformations 
which has been formulated by Green, Rivlin, and Shield (4). Stress-strain 
relations and equations of equilibrium are derived which, together with 
suitable boundary conditions and a knowledge of the state of the body B, 
make possible the determination of the additional displacements and 
stress components that are superimposed when W is perturbed. 

In section 4 the theory is applied to the problem of simultaneous ex- 
tension, inflation, and shear of a cylindrical tube. The solution of this 
problem for the case of a Mooney solid is due to Rivlin (5). In the present 
application the additional small displacements and increments to the 
components of the stress tensor that arise when the material has a strain- 
energy function of the form (1.4) are calculated in terms of the scalar 
invariants @W’/él, and @W’/éel,, and integrals involving these invariants. 
The increments to the surface tractions that are required to maintain the 
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deformation are also calculated. In section 5 the special case in which 
the tube undergoes shear but no extension or inflation is treated in a 
similar manner. 


2. Notation and formulae for finite deformations 

The notation adopted is that used by Green and Zerna (6), where 
a full account of the derivation of the theory will be found. For con- 
venience, the principal formulae and notations we shall employ are sum- 
marized here. 

The points P, of an undeformed elastic body Bo, at rest at time t = 0, 
are defined by a system of rectangular Cartesian coordinates 2;, or by a 
system of general curvilinear coordinates @;. The line element ds, of By 
is given by 


ds? — dx, dx; — 4,40; dO,.. (2.1) 
J, is the covariant metric tensor of the unstrained body By. The contra- 
variant metric tensor of B, is denoted g“ and is defined by 

Ire = Bee 
The usual summation convention is used and indices take the values 1, 2, 
and 3. The Kronecker delta is denoted by 38}.. 

The body B, is now deformed into a strained body B, so that at time ¢ 
the material particle which was originally at the point P, has moved to 
the point P. The curvilinear coordinates 6; move with the body as it is 
deformed, and form a curvilinear coordinate system in B at time ¢, in 
terms of which coordinate system the line element ds in B is given by 

ds* = G;,,d0,d6,,. (2.3) 

G,,, is the covariant metric tensor of the strained body B at time t. The 
contravariant metric tensor of B at time t is G“, where 

GG, = dj. (2.4) 

The vector P,P is the displacement vector, denoted v(6,,6,, 6,t), and 
r(@,,0,,0,) is the position vector of P, referred to the origin of the 2; 
coordinate system. The covariant and contravariant base vectors G;, G‘ 
at P of the coordinate system 6; are then defined by 

G; = r,t Viv Gi = GkG,, (2.5) 

and we have G,;.G, = Gy, G'.G* = Gi, (2.6) 

A comma denotes partial differentiation with respect to the 6, coordinates. 
The most convenient form of the equations of motion is 


T+ pF* = pf*, (2.7) 
where +“ is the contravariant stress tensor referred to the coordinates 0,, 
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p is the density of the strained body B and F* and f* are the contra- 
variant components of the body force vector and acceleration vector 
respectively, both measured per unit mass of B and referred to the base 
vectors G,. The double line denotes covariant differentiation with respect 
to B, that is, with respect to @' coordinates and the metric tensor com- 
ponents G*, G,,. The Christoffel symbols of B are given by 
‘te = 8G (Gye t+ Ges— Fires) (2.8) 
When the body B, is homogeneous and isotropic the strain-energy W, 


measured per unit volume of the unstrained body Bp, is a function of the 
strain invariants J,, J,, J;, so that 
W = WI, 1, 4,), (2.9) 
where 
I, g°G,,, I, IO", Is Gig, G= |Gyl, 


(2.10) 
The stress-strain relations then take the form 


rik — Ogik +P Bk + pGik, (2.11) 
2 0oW ow 

: eee eee — 2J/i-_., 
It r I, ’ Pp iy y 


Bik — (gikg*—girg**)G,,. (2.13) 


where = — — Y = 


If the material is also incompressible, then G = g, and J, = 1, at all 
points of the body, and W is a function of J, and J, only. The stress- 
strain relation (2.11) remains valid, but now 
2°”, y= 22 

el, el, 
and p is a scalar invariant function of the coordinates 6;. 

When the surface force vector P is prescribed at a surface with unit 
outward normal vector n, then at that surface 


® (2.14) 


rikn, — Pk, (2.15) 


where P= G,, n = n,G'. 


3. General theory 

We will now consider the effect on the stresses and displacements in 
the body B of replacing the strain-energy function W of the previous 
section by a perturbed strain-energy function W+e«W’. It is assumed 
that, with some given strain-energy function W, a solution of the equa- 
tions of section 2 has been obtained for some body, which in its initial 
and first deformed states is denoted by B, and B respectively. This is to 
say that the displacement vector v and the stress tensor 7**, which describe 
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the deformation and the state of stress of the body B, have been deter- 
mined as functions of a set of curvilinear coordinates 6; and the deriva- 
tives of W(/,, /,, J), and that these quantities, together with the metric 
tensors and base vectors of B, and B, may be regarded as known in all 
that follows. 

Now suppose that the strain-energy function W is changed to a new 
strain-energy function W-+«W’, where « is a small constant. This will 
result in a further deformation of the body B into a second deformed body 
which will be denoted B’. In this further deformation, the points which 
were at P, in By, and at P in B, are displaced to P’ in B’. The prescribed 
surface tractions on the surface of the body B’ are allowed to differ by 
quantities of at most order « from the values they have on the surface 
of the body B. It is assumed that under these conditions the incremental 


displacement vector PP’ will be of order ¢ at every point of the body, so 
that we may write 

P,P’ = P,P+ew = v(6,,t)+ew(O,,t), (3.1) 
and in the following analysis neglect «? and higher powers of « in com- 
parison with unity. 

As far as the treatment of the geometry of the deformation is con- 
cerned, the theory which follows is identical with the theory of small 
deformations superposed on finite deformations which has been described 
by Green, Rivlin, and Shield (4). We may therefore make immediate use 
of certain results of this paper, which are summarized in the following 
paragraph. 

The covariant base vectors of the coordinate system @; at points P’ of 
the body B’ are denoted by G;-+-«G;, where 

Gi = wy. (3.2) 

The displacement vector w is most conveniently expressed in components 
referred to the base vectors G,, G‘, as follows: 

w = w,, G" = w"G,,, (3.3) 

so that G; = w,,||;G" = w™|,;G,,. (3.4) 


The covariant metric tensor of B’, evaluated at time t, is G,,+-«G,, where 
Gi, = Wille Hwee, (3.5) 

and the contravariant metric tensor of B’ at time ¢t is G*+-«G’*, where 
Gk — —Girqeg,. (3.6) 


The determinant of the metric tensor components G,,-+-«G, is denoted 


by G+ eG’, where G’ = GQ*G,. (3.7) 
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Finally, the strain invariants associated with the body B’ are denoted by 
I,+el,, I,+€1, I,+<¢13, where 
T= G9"°Gre Ty = Gy G1, +G"I3), Is = G'/g = 1, G"G,,,. 
(3.8) 
For the body B, the strain-energy function has the form W(J,, J,, J;). 
For the body B’, this becomes 
W(L, +1}, I,+-e1g, [g+-e1g)+«W' (+e, 1,414, I5+-€15). 
The scalar invariants associated with the body B are ®, ’, and p, and are 
defined by equation (2.12). The scalar invariants in the body B’ will be 
denoted by D+’, Y’+-e'¥"’, p+ep’. These scalar invariants in B’ are each 
derived partly from the strain-energy function W, and partly from the 
strain-energy function W’. The terms depending on 
W (1,4 el}, I+ ely, I, +I) 


are the same functions of J,+-e/}, J,+-el3, and J,+«/, as ®, ‘Y’, and p are 
of J,, f,, and J,, and so, by Taylor’s expansion to order ¢, these terms 
are 


. ——— = sr) 


inp t feayp tsar 
el, él, él, 
, or sor) 


ay” ae 


el Po PP 
F (ia, ‘of, “*al, 





To these must be added the invariants derived from the perturbing strain- 

energy function «W’([,+-¢J}, [,+ 1, 1;+«I5). To the first order in «, 
these terms are 

2 aw’ 2 ew’ é 

€- «2/3 


“he OR 


Ww’ 


- 3.10 
aI; (3.10) 


Combining the terms of order ¢ in (3.9) and (3.10), and making use of 
(2.12), we have the final expressions for ©’, ‘Y’, and p’, viz. 

® . , 3 av’ 

9 I; } o 
2], I} él, 

. r 2 aw’ 

21, * Ty al, 

- (£1, +dr,+cr,) +. r4+2n2@ 

2/, él, 


®’ = Al, + FI,+£EI,— 
FI, + Bly+DI,— (3.11) 


) 
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2 @w 2 ew 2 aw 
ay2? B=—-; 2’ ~ Ti er 


ew 2 &w 2 &w 


A —> 
(3.12) 
D = 


" al,ai, ~ Noah,’ ~~ Rahal, 


For an incompressible body, J, = 1 and J; = 0. W and W’ are func- 
tions of J, and J, only, so that 
OW 


= 257 (3.13) 
and p is a sealar function of the coordinates which is determined by the 
equations of motion and boundary conditions for the body B. Also A, 
F, and B may still be found from (3.12) by setting J, equal to unity, and 
we have re 
©’ = AI,+ Pr,+ 255 
pads 3 (3.14) 
W = FL,+BM,+2o 
él, 
The pressure function p’ cannot now be derived from the elastic potential 
and is a scalar invariant function which must be determined by the 
equations of motion and boundary conditions for the body B’. 

The remainder of the theory is again formally identical to the theory 
of small deformations superposed on finite deformations, but it should be 
noted that the functions ®’, ‘Y’, and p’ of the present theory differ from 
the functions ®’, ‘’’, and p’ which appear in the paper by Green, Rivlin, 
and Shield by the addition of terms involving the perturbing strain-energy 
function W’. 

The tensor B becomes, in the body B’, a tensor B*-+-«B’, where 


Bk — (g'*g"*§—g'"g**)G,. (3.15) 
To the first order in ¢ the stress tensor for the body B’, referred to 


curvilinear coordinates 6; in B, is r-+-er'*. By substitution in the stress- 
strain relations, we find 


7th = gtk’ +. Beep’ +. Ben + Gtk’ 4 G’iky, (3.16) 
If the contravariant components of the body force and acceleration 
vectors for B’ are respectively F*+-<F’* and f*+<«f", referred to base 


vectors G,+-«G;,, and p is the density of B, then the equations of motion 
for B’ are 


[7 rtkwr ||, +r ||) |p pL P'* + Fret ||,] = pl f’*+frw*||,]. (3.17) 
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In the case in which surface forces are prescribed on B’, and n;,+en; 
are the covariant components of the unit normal to the surface of B’, 
referred to base vectors G'-+-«G"', while P*+«P” are the components of 
the surface force vector, then the boundary condition at the surface of 
B’ is ; 

(r*+er’t*)(n; +en;) = Pk+eP"*. (3.18) 


4. Shear, extension, and inflation of a cylindrical tube 

As an example of the application of the preceding theory, we now 
consider the simultaneous shear, extension, and inflation of an incom- 
pressible cylindrical tube. This problem has been solved in the case in 
which the strain-energy function of the body has the Mooney form by 
Rivlin (5), and by a rather different method by Green and Zerna (6). 
We use here the method and notation of Green and Zerna. 

A circular cylindrical tube which in its unstressed state has external 
and internal radii a, and a,, and length /, is subjected to the following 
successive deformations: 

(i) A uniform simple extension of extension ratio A. 

(ii) A uniform inflation of the tube, in which the outer and inner radii 
of the tube change to r, = p, a, and r, = ped,, respectively, the length 
remaining unaltered. 

(iii) A shear of the tube about its axis, in which each element of the 
tube rotates about the axis through an angle ¢ which depends only on the 
radial position of the element. 

(iv) A shear in the direction parallel to the axis, in which each element 
moves in a direction parallel to the axis through a distance w which 
depends only on the radial position of the element. 

The moving coordinates (@,,4,,@,) are chosen so that they coincide with 
polar coordinates in the deformed body. The point (r, @,z) in the deformed 
body was originally at the point (aw, , ¢), where 


w=-aolr)—rQr), @=g4+d(r), z= dAlC+wir), (4.1) 
and hence, in the notation defined in section 2, 


2 rQ(r)cos(@—¢), x, = rQ(r)sin(@—d), x3 = [z—wi(r)]/A. (4.2) 


The theory outlined in section 2 may be applied without making any 
assumption regarding the form of the strain-energy function, as far as 
the formulation of the equations of equilibrium of the deformed cylinder. 
This is done by Rivlin and by Green and Zerna. The following results 
are quoted from Green and Zerna. The metric tensors of the undeformed 
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body are r 2 


Gi +rQrg? +5 
—rQd, 
= 
d2 
® a "6, 

A? ? 
Od, 1 | Vor 
2 ry? v 2 
ky Q*4, wv, 
A? A? 














where the suffix r denotes differentiation with respect to the r coordinate. 
The metric tensors of the deformed body are 


1 0 0 . ° 
Gy = |0 rr O}, G*=10 0 


0 0 1 


(4.5) 
0 1 


In equations (4.3) and (4.4) the expressions for the components of g,, and 
g‘* have been simplified by use of the incompressibility condition 
g=G=r, 
which determines Q as a function of r by the equation 
rQ = {APr+K)}, 
where = a3(1—Au?) = a}(1—Ay§). 


The strain invariants are 


20242 2 
h=Gt¥+h ee 


I, = +54 Mtoe st 


and the tensor B* is 

1 

2 + Q? "4, 
an 


Q%, scat jyat OMt+ 


— 


+58 


OMT 








| 0 P+a 


Q? 
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The components of the stress tensor are then shown to be 


2 ] r, TORN 
oo (Lie)rip 


P ] r2()*4?) (= l —.; 
= » abe r = 1m x 3 j =( 22 1 r y Lay 
(oi » je qty t eer ") i 


= (184 Selo 4 (45 + eggs SAN +p 
ry Q ® 
; (4.10) 


= Cre 4.044, 


Q?w Wve 
a rp ry 


7 


7 





723 Q¢, vr® 
2 


The non-zero Christoffel symbols are 
aS es 
ri. = —r, ri, = Th, -. (4.11) 
r 
Assuming zero body forces, and noting that the components of the stress 
tensor are functions of r only, the equations of equilibrium are found 


ell 
dz % ] (ql p2p22y — ¢ 
dr 


to be 





The last two of these may be integrated immediately, and on substituting 
for 7! and 7!% they yield 


2B, 
~ BQO +P) | 


2D, 
r(Q°o+¥) 


where B, and D, are constants of integration. 

In general, ® and ¥’ are functions of J, and J,, which are shown by 
equation (4.8) to be dependent on ¢, w, and r, so that equations (4.13) 
are non-linear differential equations which determine ¢ and w as func- 
tions of r. To integrate these equations we must assume some definite 
form for the strain-energy function W. If we suppose W has the Mooney 
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form, defined by equation (1.1), then 
® = 2, Y = 2,, 
and the equations of equilibrium can be integrated to give 
feat AK(C. eer 5+ +8, 


D, 5 (CG, QP+€,)r°) 
oinlei I sat PR el 
AC, 0) 8G, QeLC,yra| + 2 
Bi | 
MK>AC,+2°C,)| 
D30, 
~ NrHC, G+-C,)AC, +C,) 





p mengemys 








20. 
ree 


where Q,r, = a, F is an arbitrary constant hydrostatic pressure, and 
B,, B,, D,, and D, are constants of integration which are determined by 
the conditions at the curved surfaces of the tube. For example, if 
w=d=O0atr=—r,, and w = wp, d = dat r =r, we find 


p, — >K(G+¥C)bo ap _ 900,40, » oe (C +0, vai} 


log(Q1/Qz) \(C, Q3+-C,)r3’ 


B, = D, = 0, (4.15) 
where Q,7r, = 4. 

The surface tractions may now readily be determined by the applica- 
tion of equation (2.15). The expressions for the surface tractions are 
given by Rivlin and by Green and Zerna. 

Now consider the modifications that must be made if the strain-energy 
function W is changed to W+«W’. An additional deformation of the 
body will then in general be superimposed on the previous deformation. 
As the radial symmetry of the body is unaltered, this additional deforma- 
tion will depend only on the coordinate r. We suppose that in addition 
to the previous finite deformations (i), (ii), (iii), and (iv), the cylinder now 
undergoes the following small deformations: 

(v) A shear of the tube about its axis, in which each element rotates 
through an angle ed’, which is a function of ronly. ~* 

(vi) A shear parallel to the axis, in which each element moves through 
a distance ew’, which is a function only of the radial position of the 
element. 

Since Q(r) is shown by (4.6) and (4.7) to be independent of the form of 
the strain-energy function, perturbing the strain-energy function causes 
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no additional displacements in the radial direction provided that the 
internal and external radii of the tube are kept constant. The displace- 
ment vector w therefore has components in the @ direction, determined 
by (v), and in the z-direction, determined by (vi), but none in the radial 
direction. 

With the above assumptions about the form of the deformation, the 
displacement vector w has, in the notation of section 3, the following 
components referred to the base vectors of the polar coordinates (r, 6, z) 
in the first deformed body: 


w! = 0, w? = ¢'(r), w> = w'(r) 


w, = 0, WwW, = r*¢'(r) 
Hence, from (3.5) and (3.6), 


0 rd, w, 


1: = rg) 0 0 
w. OO O 


and from (3.8) 


f= Os 6 See 


21 Qu*4,4,+-5 > 3 Yr w | 


0 


where the suffix r again denotes differentiation with respect to r. It 
follows, from (3.14) that 


c 


®’ = 2A Se “Tt d, .@ ~ = W, w, | +2F | a4, d+ a w, u,| +2 wil \ 
(4.19) 


2 2 oy’ 
2F| SS d, >, @ w, 1, +28 724 db, +7% ui] +257 
r? ? : él, 


Substituting in (3.15), the tensor B’* is seen to be 


$, : 
=o 2 sisi Qu r 


’ 
W, WwW, 


1 , , 
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and, applying (3.16), the incremental stress tensor 7’ has components 


= Q* 5 p ) 
ll 2 
= 20+ (a+ @e)e +p 
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Since by symmetry the stress components are functions of r alone, the 
equations of equilibrium (3.17) become 

: [rt ppteg’] 4 | ppt ptp’22_ 919 Ong’ 4 8g’) =0 
dr r 


d 12 ufyr ty, 1 2 only, tar) peer] _ 
5 [rtrd +o) | +2 [arte +20n(4, + 24°) reg] = 0 
d [7's + 2M] + =f’ + ates] =? 


dr 





Making use of the first two of equations (4.12), the equations of equi- 
librium (4.22) may be rewritten 
dr’ ) 


oo + Lf tet] arettgy = 0 


3 
© 721121 wlld’ [7124-711 6’) — 
S [rt agi] 4 Spe 24 og’) = 0 
d 
= 
The last two of these may be integrated immediately, giving 
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where B;, and D;, are constants of integration. 
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At this stage it is necessary to assume a particular form for the strain- 
energy function W. We may, however, proceed without making any 
assumption about «W’, other than that ¢ is small. If we suppose that W 
has the Mooney form, then 


® = 2C,, "= 20,, 


- m 
o’ — 2° Ww’ yr — 2° Ww’ 
& b ad > 
el, él, 
and ¢, w, and p are given by equations (4.14). Equations (4.24) then 
become, after substituting for the stress components, 


, l 2 a oy LAMe” | 
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Te 
where Bi, and Dj are constants of integration. 
If we require as boundary conditions that the displacements at the 
curved surface of the second deformed body are to be the same as the 
displacements at the curved surface of the first deformed body, then 


¢ =w' =0 when =7r, andr=r,, 
and it follows that EK, = D, = 0, 
QO 


B q dQ 
a 4 sure 
. 2A7(C, +A*C, log (Q,/ Qs) | +A") Q 
OQ: 


Di, = | 2D,(C,+2°C, ) fr? og! | C at f (Q' -¥") dr 
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The first of equations (4.23) then gives, after integration and some 
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simplification, 
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\¢+ F’, (4.28) 


where F’ is a constant which represents an arbitrary hydrostatic pressure. 
If the equation of the surface of the body B’ is 


F (4,95, 5) = 0, (4.29) 


then the covariant components of the unit outward normal vector n to B’, 
referred to base vectors G'+-«G", are 


4 


n, ten, = ka 
where & is chosen so that n is a unit vector. In the present example, the 
curved surfaces r = r,, r =r, in B become, in B’, with the boundary 
conditions that have been adopted, the surfaces #, = r, and 6, = r,, so 
that n has components (1, 0,0) at the outer curved surface of the cylinder, 
and (—1,0,0) at the inner curved surface. Equation (3.18) then gives as 
the components of the surface force vectors 


(4.30) 


Pk+A- PRP" (r*), nte(r*), om at 6, == r:) 
, (4.31) 
PkA-¢ P’* —(r"*), net’) par, at 6, =" = ry) 


these components being referred to the base vectors G,;+-«G;. 
The ends of the cylinder, which in the body B, are plane, are deformed 
in the body B to the surfaces 


z+Al—wir) = 0 
and in the body B’ to the surfaces 
z+Al—w(r)—ew'(r) = 0. 

Since the coordinates 6; in the body B’ are related to the coordinates r, 0, z 
by the equations r= 06, 

6 = 6,4<«¢’ }, (4.32) 

z= 0,+ew" 
the equations of the end surfaces of the cylinder may be written in the 


form 6,+A—w(8,) = 0. (4.33) 
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On the surface which corresponds to the positive sign in equation (4.33), 
the vector n has components 


(1+-w?)!(n,-+€n}) 


, 


(1+-w?)#(n.+- eng) 
(1-+-w?)t(n,+en3) = 
and so at this surface 
(1-+-w?)?( P*+¢P*) = —w,(rl*+ er’) + (73*+ 7%), (4.35) 
where again the components of the surface force vector are referred to 
base vectors G,;+-«G;,. 

Alternatively we may express the surface force vectors in components 
Q*+.«Q’* referred to base vectors G;. In order to do this we transform 
from the coordinates (@,,0,,@,) to the coordinates (r,@,z), using the rela- 
tions (4.32). On making this transformation, we obtain 

Q+«eQ" = Pi+eP" 
Q?4 «Q” ” P?+-«(P’*+¢4) P') . (4.36) 
Q?+€Q’ = P3+¢(P%+w;, P?) 


5. Shear of a cylindrical tube 


We now consider the special case of the problem of the preceding 
section in which the tube is not extended or inflated, so that 


Hy = Hg =A l. 


It follows that @ = Q, = QY, = 1, and K = 0. In this case, again follow- 
ing Green and Zerna, the solution of equations (4.13) applied to a Mooney 
solid, with boundary conditions ¢ = w = Oatr = a,, and ¢ = dy, w = wy 
at r = a; is 


#1 D, log(r/az) 
C,+0, 


w 


(5.1) 


and the pressure function p is found to be 


B? C, D? 
2(C,+C,)r* (C,+C,)?r? 





p= +F, 


ct 2(C,4 C,)aj a3 bo - s (C,+C,)w, 


where B, 





-C,)a = Cat la)ieo 5.3 
a?—az . log(a,/a.) "o 


The incremental stress components 7’ are again given by (4.21), with 
appropriate substitutions for the special case now being considered, and 
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as before equations (4.23) are the equilibrium equations for the second 
deformed body. Equations (4.25) now simplify to 
; l 0’ +" 
eee a 
(( =o © *G+G, 
l 0’ +4" 
nnn | PY — 3) ——__ 
orl ve ater 


which on integration give 
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a2 
If we again impose the conditions w’ = ¢’ = 0 when r = a,, r = ay, then 
B= D,=0 ) 


ay, 


aiaz B, [ o's ea 
_ SS fea 
@i—anc,ta) J Ott 
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BY 


[oiw)% 
: r 

Ie } 
p’ is then found from the first of equations (4.23) to have the value 
B, B, 2D, DC, 


iC, +0.) FC,+6,7 


dD, 


Dy : a 
2 log(a,/a,)(C, +C,) 





p' = F’—@’—2¥" 





"| BE @ +8) Dt [20,@'—(C.—G¥"Jdr 
\(e,4+¢,? rf "(4067 r2 ro 


(5.7) 


Equations (4.31) and (4.35) again give the components of the surface force 
vector. 
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EXACT SOLUTIONS FOR THE GROWTH OF FINGERS 
FROM A FLAT INTERFACE BETWEEN TWO FLUIDS 
IN A POROUS MEDIUM OR HELE-SHAW CELL 


By P. G. SAFFMAN (Cavendish Laboratory, Cambridge) 
[Received 6 February 1958] 


SUMMARY 
A family of exact solutions of the equations of motion for the two-dimensional 
unsteady motion of a viscous liquid driven by a fluid of negligible viscosity in a 
porous medium or Hele-Shaw cell is obtained. The solutions represent the growth 
of a particular disturbance of a state of uniform motion and its development into 
a motion in which long fingers of arbitrary width and spacing penetrate the viscous 
liquid. 


1. Introduction 

IN a recent paper, Saffman and Taylor (1) have considered the motion of 
the interface between two viscous fluids in a porous medium or Hele-Shaw 
cell when one of the fluids is driven by the other. The flow in a Hele- 
Shaw cell consisting of two parallel, closely spaced, flat plates is mathe- 
matically analogous to two-dimensional flow in a porous medium in which 
the motion is governed by Darcy’s law. 

The stability of a steady state of uniform motion, in which the interface 
is plane and perpendicular to the direction of flow, was investigated both 
theoretically and experimentally and it was shown that small disturbances 
of the interface grow if the less viscous fluid drives the more viscous with 


a velocity greater than a certain value, The later stages of this instability 


were investigated using a Hele-Shaw cell and were found to consist of 
the penetration of ‘fingers’ of the less viscous fluid into the more viscous 
one. A natural idealization of this phenomenon is to consider the propaga- 
tion of equal and equally spaced fingers, and since the streamline midway 
between two such fingers is straight and may be replaced by a channel 
wall (at least for the purposes of mathematical analysis), this led to a study 
of the steady propagation of a single long finger of fluid through a parallel- 
sided channel containing more viscous fluid. The latter motion was studied 
in some detail in (1). 

During the course of this work, exact unsteady solutions of the equations 
of motion were discovered for the case in which interfacial stress effects 
are negligible and the viscosity of the driving fluid can be neglected. The 
solutions represent the growth of long fingers propagating steadily through 
the viscous liquid from a flat interface which is given a particular initial 
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disturbance. The form of the initial disturbance is approximately sinusoidal 
when its amplitude is small. The motion is periodic with arbitrary wave- 
length in a direction perpendicular to that of the velocity at infinity, and 
the flow in one complete wavelength represents the growth of a single 
finger in a parallel-sided channel. The solution does not appear in itself 
to be of much physical importance, since the particular initial conditions 
cannot be exactly reproduced in practice. However, it is of some mathe- 
matical interest, being an exact solution of a non-linear problem, and does 
have some bearing on a problem of physical interest. 


2. The equations of motion 


We consider here the motion of a viscous liquid which is driven by a 
fluid of negligible viscosity in a region of infinite extent. The theory is 
restricted to the case in which the two fluids do not interpenetrate to any 
marked extent and can be regarded as separated by a sharp interface; 
and the analysis applies to two-dimensional motion in a porous medium 
or flow in a Hele-Shaw cell. For definiteness we consider the latter. Let 
x and y be coordinates in the plane of the plates bounding the cell, with 
the x-axis taken parallel to the velocity at infinity; then neglecting the 
effects of gravity (these can be taken into account without difficulty but 
they complicate the algebra and it is found that they do not alter the 


character of the solutions in any way), the mean velocity across the 
stratum of viscous fluid contained between the plates is given by 


b? op op  — oo 


— , 


A 


l2u éx Ox cy 
onsite ee (1) 
l2u Cy) Gy x 

where u and v denote the components of mean velocity parallel to the 
x- and y-axes respectively, p is the pressure, » the viscosity of the liquid, 
b the gap between the plates, ¢ = —(b?/12u)p is the velocity potential, 
and the stream function ¢% exists by virtue of the equation of continuity. 
The viscous liquid is taken as extending to x = +00 and the other fluid 
to x = —«. The velocity at infinity is supposed uniform and is denoted 
by V; hence ¢ ~ Vz as x > +00. 

The pressure is constant in the fluid of negligible viscosity, and if effects 
due to interfacial stress can be neglected (see (1) for a discussion of this 
assumption), the pressure in the viscous liquid is constant along the inter- 
face. The equation of the interface can therefore be taken as ¢(z, y,t) = 0. 
We assume further that the viscous liquid is completely expelled by the 
other (this assumption is also discussed in (1)), so that the interface is 
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a material surface moving with the mean velocity; hence 


Op , 9h, % 


Cc C2 cy 


0 


on ¢ = 0. 

The position of the interface in the physical plane is not known a priori, 
but it is specified in the potential (¢,%) plane. We therefore work in the 
potential plane and consider x and y as functions of ¢ and y. The region 
occupied by the viscous liquid corresponds with the domain ¢ > 0. 
Further, ¢ and & satisfy the Cauchy—Riemann equations, and therefore 
z = x+y is an analytic function of w = 4+, satisfying 


9 V as db > +7 (3) 


because the velocity at infinity is uniform. The boundary condition on z 
at the interface is obtained by expressing (2) with ¢ and & as independent 
variables. A straightforward application of the rules for change of variable 
in partial differentiation gives 


Ox Cy Ox ey 
et xs xb Ot (#) 
on @ 0. 

Any function z(w), which is analytic for 6 > 0 and satisfies (3) and (4), 
gives a mathematically possible flow. (4) is non-linear and the problem 
of finding a solution to fit given initial conditions is in general difficult. 


3. Exact solutions representing the growth of fingers 


The expression 


7 | me. 


where V = AU and —V/ < & < VI, was obtained in (1) for the steady 
propagation with velocity U of a semi-infinite finger of asymptotic width 
2\/ in a channel bounded by walls y +l containing viscous liquid 
moving with velocity V at infinity.t If % ranges from —oo to +, (5) 
represents the motion of equal and equally spaced semi-infinite fingers, 
each of width 2A/ and separated from its neighbour by a gap of width 
2(1—A)l, propagating steadily through the viscous liquid. It is clear (as 
pointed out in (1)) that there exists a family of solutions corresponding 
with values of A between 0 and 1, and the case of steady motion does not 
possess a mathematically unique solution. 


= (1—Allog 4{1 +-exp/ vi) (5) 


Whilst searching for a mechanism to distinguish a particular value of 


+ The term Ut, which is not present in the equivalent expression given in (1), arises here 
because the origin is fixed in space, whereas in (1) it was fixed relative to the finger. 








EXACT SOLUTIONS FOR THE GROWTH OF FINGERS 
A and make the solution physically unique, it was discovered that 


z = (+ dtl) += (1—djlog 41 +a(texp( —F7)} (6) 
J 7 \ vi}| 


is analytic ind > 0if 0 <A < land —1 < a(t) < 1, and satisfies (3) and 
(4) if (1—A)lé+-mdad = 2Va, 
2(1 —A)(2A—1 lad +-n{1+-(2A—1)a*}d = 2V(1+-2%), (7) 


where the dot denotes differentiation with respect to time. The integrals 
of (7) are easily found to be 
a 
(1- a2)2Aa A) ~ 


yerV Ul d = Vt— 2) _))*log(1 —a*)+8, (8) 
7 


where « and £ are arbitrary real constants. 
The interface between the two fluids is the image of ¢ = 0, and its 
parametric equation is 


x = d(t) “(1—Ajlog if +.a# +24 000( 7) 


a Vij)’ 


yan af _@sin@ed/VD_ 
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y (9) 


7 
Examining the shape of the interface as given by (9), it is to be noticed 
in the first place that there exists a whole family of possible shapes, 
corresponding with values of A between 0 and 1, all of which are periodic 
in y with period 2/. If a< 1 and t < (l/V )log(1/«), then the interface has 


the equation a = Vt+A cos(ry/le"™4+ O(a2), (10) 


where A = (2l/m)a(1—A). This represents a sinusoidal perturbation with 
amplitude A of the flat interface in the state of uniform motion, which 
grows exponentially with amplification factor 7Vt/l. It is to be noticed 
that the form of the interface is, to the first order in the amplitude, 
independent of A. Further, as t > 00 


1—a ~ exp[ —aV1t/{2/A(1—A)}], d ~ Vt/A, (11) 


from which it follows that the motion given by (6) tends to that given by 
(5). Thus, for large ¢ the interface consists of equal and equally spaced 
long fingers, of width 2A/ and length of order Vt/A, propagating into the 
viscous fluid, the tips of the fingers advancing with velocity V/A. The 
solution thus represents the growth of long fingers of arbitrary width from 
particular, approximately sinusoidal, disturbances of a flat interface. 

In Fig. 1 are shown one complete wavelength of the interface for three 
pairs of values of A and a. Curve (a) corresponds with A = 0-8, a = 0-95; 
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curve (6) with A = 0-2,a = 0-5. The amplitudes of the two curves are 
roughly the same, about } of the wavelength, but the flatness of that 
corresponding to A = 0-8 is already pronounced. Curve (c) shows the 
interface with A = 0-5 when the amplitude is somewhat larger; for this 
curve a = 0-95. The dotted line shows the position of the interface when 
there is no disturbance. 








Fic. 1. Shapes of one complete wavelength of the interface for different values 

of A. (a) A 0-8, a 0-95; (b) A 0-2, a 0-5; (ec) A 0-5, a 0-95. 

It was found experimentally in (1) that the widths of fingers grown from 
a flat interface were of the same order as the distance between them. 
Experiments with single fingers in a channel showed that A was close to 
} when variations in the pressure drop across the interface due to surface 
tension are apparently negligible. Now from (11), a > 1 least rapidly or, 
in other words, the exponential approach to the final asymptotic state is 
slowest, when A = }. This distinguishes a particular value of A but the 
physical significance, if any, of this result is not clear. No other feature 
of the solutions which throws light on the occurrence of A = 3 in practice 
has been detected. 
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THE AERODYNAMIC FORCES ON AN AEROFOIL IN 
UNSTEADY MOTION BETWEEN POROUS WALLS 
By 8. ROSENBLAT 


(Department of Mathematics, University of Manchester) 


[Received 2 January 1958] 


SUMMARY 

A solution is obtained to the boundary-value problem arising in the unsteady 
motion of a thin aerofoil in a stream bounded by porous walls. The flow is two- 
dimensional, inviscid, and incompressible. The condition holding at the porous 
boundaries is assumed to be a proportionality between the normal velocity com- 
ponent at the wall and the pressure difference across the wall. The solution relies 
on a form of this boundary condition which is valid for small values of the frequency 
of the oscillations, although conditions at the aerofoil surface are independent of 
frequency. The usual assumption of small amplitude of the disturbances is made. 

A formula is derived for the pressure on the aerofoil surface in terms of the flow 
direction at the surface. The relations for the lift and moment about the mid-chord 
point are obtained for a harmonic upwash distribution. For the case of a rigid 
aerofoil, particular types of harmonic oscillations are investigated, and the lift and 
moment are given as dimensionless ‘air-load coefficients’. The results are expressed 
as the first-order terms of expansions in powers of a parameter which depends on 
the ratio of aerofoil chord to tunnel width, and are valid for small values of this 
ratio. 

1. Introduction 

Tue problem of an aerofoil moving in a tunnel with porous walls has 
received comparatively little attention. In a recent paper (1), Drake 
obtained an extended form of Possio’s integral equation for an oscillating 
aerofoil in two-dimensional compressible flow between porous walls. For 
the steady flow problem, Baldwin, Turner, and Knechtel (2) have pre- 
sented an approximate solution for lift and drag interference in two- 
dimensional and axially-symmetrical tunnels with both porous and slotted 
walls. Further solutions for drag interference and the steady lifting aero- 
foil have been given by Woods (3), (4). 

Particular cases of a porous-wall tunnel are those in which the porosity 
is zero (solid-wall tunnel) and infinite (free jet). Solutions to these problems 
have been given, for the unsteady aerofoil, by Rosenblat (5) and Woods 
(6) respectively. The present paper is a generalization of these works, 
and the analysis follows in the main the pattern set by them. 

The basic difference between the present problem and the solid-wall 
and free jet problems lies in the nature of the boundary condition at the 
walls. For the solid-wall tunnel the component of velocity normal to the 
walls is constant, while for a free jet the pressure is constant along the 
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jet boundaries. In the present case it is assumed that the velocity com- 
ponent normal to the wall is proportional to the pressure difference across 
the wall, where the constant of proportionality is a measure of the porosity. 


This relationship has been used in their work on porous-wall tunnels by 
the above-mentioned authors. 

It is shown in section 2 below that for oscillations of small amplitude 
and frequency this condition at the boundary can be expressed as 

0 = wQ, 
where @ is the flow direction at the walls, w is the ‘porosity’ of the walls, 
Q is the ‘velocity logarithm’ = In(U’/q), q being the fluid velocity and U 
a standard reference velocity. 

The aerofoil is taken to be situated midway between the walls, and 
the usual assumptions of unsteady aerofoil theory are made. That is, we 
suppose (i) that the aerofoil is sufficiently thin for the two components of 
the flow direction at the surface, that is, the symmetrical part associated 
with the unsteady upwash and the anti-symmetrical part arising from 
thickness effects, to be additive—the latter part being neglected in the 
integrations as it makes no contribution to the lift and moment; (ii) that 
an infinitely thin ‘wake’ or vortex sheet is generated from the trailing 
edge; (iii) that the displacements due to unsteadiness are small enough to 
allow only first-order terms in the perturbation velocity to be taken. 

The physical plane of the problem is conformally mapped, by a trans- 
formation involving Jacobian elliptic functions, into a rectangle of which 
one side represents the aerofoil surface, the opposite side the tunnel walls, 
and the other pair of sides the upper and lower surfaces of the vortex 
sheet. It is then seen that the function 

f =Q+70 = In(U/q)+-100 
is analytic everywhere within and on this rectangle, except possibly for 
a finite number of logarithmic singularities on its border. Moreover f is 
determined from the boundary conditions (i) known @ at the aerofoil 
surface, (ii) proportionality of @ and 2 at the walls, and (iii) continuity of 
pressure across the vortex sheet. An appropriate expression for / in terms 
of its known boundary values is found in the form of an integral equation. 

The requirement that f vanishes at infinity upstream yields an equation 
relating 6*, the flow direction at the aerofoil surface, and X, the jump in 
the velocity logarithm across the wake. A further integral condition con- 
necting these two functions is found from considering the circulation about 
a closed contour which includes within it the aerofoil and wake. The func- 
tion X is also seen to satisfy a partial differential equation derived from 
the fact that the pressure must be continuous across the vortex sheet. 
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These relations are subsequently used in the determination of the lift and 
moment. 

An expression for the pressure distribution on the aerofoil surface is 
next derived. By using the above-mentioned relations the pressure is 
given in a form not involving the function X explicitly. From the pressure, 
formulae for the lift and moment about the mid-chord point are obtained. 
Some of the integrals encountered in this determination cannot be 
evaluated exactly, and it becomes necessary to resort to approximations. 
It is assumed at this stage that the ratio of aerofoil length to tunnel width 
is small, and expansions are carried out in terms of a parameter k, which 
is a direct function of this ratio. Odd powers of k are found to vanish, 
and formulae for the lift and moment are obtained which are correct to 
order k?. 

These formulae are in terms of 6*, the flow direction on the aerofoil 
surface, which, however, is not completely specified a priori. It is shown 
that 6* is compounded of a term proportional to the prescribed upwash 
velocity, together with a term associated with the movement of the front 
stagnation point. This latter component must be determined separately, 
from the conditions of the flow. 

As a particular case we then consider a rigid body aerofoil performing 
simultaneous translational and rotational (pitching) oscillations. For this 
case the lift and moment are expressed as real dimensionless numbers, 
called air-load coefficients. From these formulae it is possible to evaluate 
the corrections necessitated by the interference due to the presence of the 
walls. The author proposes in a future publication to discuss in detail 
these interference effects for various values of the porosity and of the 
ratio of aerofoil length to tunnel width. 

In the course of the analysis we encounter several integrals involving 
Jacobian elliptic functions. For the most part no attempt is made here 
to give the method of evaluating these integrals, and the reader is referred 
to (5), where identical or similar integrals are discussed at some length. 


2. The boundary conditions 

Consider the unsteady motion of a thin, two-dimensional aerofoil 
situated midway between straight, parallel, porous walls, as shown in 
Fig. 1. Let (¢,@) be the velocity vector in polar coordinates of an incom- 
pressible flow past the aerofoil. The origin of the z (= 2-+-iy) plane is taken 
to be at the mid-chord point of the aerofoil, whose chord length is 4a. 
The width of the tunnel is 2h, so that the tunnel walls are the lines y = +A. 

It is assumed that the aerofoil is thin enough, and its unsteady perturba- 
tions about its mean position are small enough, to allow us to impose the 
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boundary conditions prevailing on the aerofoil surface over the strip 
-~2a <x < 2a, y = 0 without significant error. It is consequently further 
assumed that the (infinitely thin) vortex sheet or wake generated from the 


trailing edge will lie on y = 0, 2a<a4<a= 


AY 


C 








y H 8B 

J iG |B 
—— 

\D 








z plane 
Fic. | 


Suppose now that the pressure outside the porous surfaces is maintained 
at some constant value, p,, say. Let v, p denote respectively the com- 
ponent of velocity normal to the walls and the fluid pressure at the walls. 
Then p—p, represents the drop in pressure across the porous surface and 
the basic assumption of the theory is that the normal velocity component 
is proportional to this pressure difference. We write 


(1) 


where p is the density and U the velocity at infinity upstream. The 
quantity @ is a dimensionless parameter, referred to as the ‘porosity’, 
whose value is determined experimentally from the character of the sur- 
face. It may range between 0 (zero porosity—solid wall) and 0 (infinite 
porosity——free jet). The use of the assumption (1) is made in (1), (3) and 
elsewhere. 

At a point on the wall the flow direction is given by 


: v\ 
@ = sin (‘), 
q 
or, to first order in the perturbation velocity q—U, 


»” 


6 = 


(P—Ps) 


Hence (1) becomes (=a 
pU? 
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Differentiation, and use of Bernoulli’s equation, gives 


U2e6 =lLép 
mt ee -|3 +; <u}. 
If now we define a ‘velocity logarithm’ 


= In(U/q), 
then to first order we find 
oq —— U a 
ot ot 
06 
so that we now have — = & peed F (3) 
és 


’ 


In order to render the problem tractable by the present method it is 
1a _ eQ 
necessary to restrict consideration to the case where Vw < —. For 
0 cs 
harmonic oscillations we may write 


Q = Q,(s)e™, 6 = 8,(s)e™, 


so that (3) becomes 
eit Po __ ett ay) — a do i +o. 
és és U 


Hence the solution will be confined to small values of the frequency of 
oscillation. The boundary condition at the wall may thus be written 


6 = wf, (4) 


where the pressure p, is selected such that the integration’ constant 
vanishes. 


3. Solution of Laplace’s equation 

Suppose F(t) is a function of t (= y+») which is regular everywhere 
within and on the rectangle —2K < y < 2K,0 < n < K’ with the possible 
exception of a finite number of logarithmic singularities on the border of 
the rectangle. Suppose further that the following data are given concerning 
the character of F(t): 

(i) the imaginary part of F(t) on 7 = 0, —2K < 2K is known, 

and is equal to im /j, say; 

(ii) the imaginary part of F(t) on » = K’, — < y < 2K is known, 

im Fy, say; 

(iii) for the sides y = +2K,0<7< KR’, 

im Fx = im F_,x¢ 

and {re F,,—re F_,x<} is known. 

An expression for F(t) satisfying these conditions has been obtained 
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(7), and is given in terms of theta functions. A more convenient form of 
this solution, involving Jacobian elliptic and zeta functions, and derived 


from it (see Appendix 1) is 


F(t) 
‘n(y* —t)dn(y* —¢) 
sn(y* —?) 


lim Fal: Ly* | 


ten(y*—t-+iA')dn(y*—t-+ik’) , ‘nies 
/ / Z(y* tt-ikh ) 
sn(y* f ik‘) 
tydn(in* —t)—1 


sn(in*—t) 


* +f) 


where 
2h 
a real constant, and A and A’ are the real and imaginary quarter-periods 


of the elliptic functions. 





”) 
lr 














The physical or z-plane is mapped into a rectangle in the t = y+in 
plane, shown in Fig. 2, by means of the conformal transformation 


en(t, k) a sinh 7 (6) 
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where k and k’ are the moduli of the elliptic functions and are given by 


7a 


k — tanh- 
h 
(7) 
or, since k?+-k’? . sinh ~~ 
/-# %e K*-+-K * : § 
k’ h 
In the t-plane the aerofoil surface y = 0, —2a < x < 2a maps into the 
line » = 0, —2K < y < 2K, the upper and lower surfaces of the vortex 
sheet map into y +2K,0< 7» < K’ and y = 2K, 0 <7 < K’ re- 
spectively, while the upper and lower walls become » = K’, 0< y< 2K 
and » = K’, —2K < y < Orespectively. The point upstream at infinity, 
F, becomes t = ik’, while the point downstream at infinity, 2, becomes 
t t+ 2K+iK’. 
If now we define a function f(t) by 
f = In(U/q)+00 = Q+7, (8) 
then f will be an analytic function of both z and ¢, except at possible 
logarithmic singularities corresponding to stagnation points on the aero- 
foil surface. Moreover it is clear that f(t) satisfies conditions (i) and (iii) 
above. For on » = 0, —2K < y < 2K the imaginary part of f is the 
prescribed flow direction at the aerofoil surface, #*, say; while for the 
sides y -2K,0< » < K’, we have that the flow direction is the same 
on either side of the wake, and that the jump in Q, which we write as 
X =Q,,—Q_.x, (9) 
is determined from the condition of pressure continuity across the vortex 
sheet. 
The conditions on f(t) may be summarized as follows: 
On » = Oand y = K’, —2K <y < 2K 
O—a = Wy) 
where on 
2K o= 0 and yb = 6* 
0 o—--—omw and #=0 
< 2K v=o and ¢=0 | 





On =O. |}. (11) 
0.56. —Q_ 57 = X 
It follows from equation (10) that a solution for f(t) cannot be obtained 


by direct application of the formula (5), since im f,. is not known. How- 
ever f(t) is derived by the following method which utilizes equation (5). 
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Let f,, = Q,,+70,, be an analytic function of t, regular in 
2K : y : 2K, 0l nS K’, 
which satisfies the ‘homogeneous’ boundary conditions: 


On 7 = Oand y = K’, —2K <y < 2K 


where on 


0, 2K <y< 2K and @ 


m 


K’, —2K <y: - and @ 


m 





O<y< 2K and @,, 


”) 
6,,(2K +n) = 6,,(—2K+in) = 0 \ 
Q,,(—2K+in) = Q1, — 
We then show that the function defined by 
® = fifm (14) 


satisfies all the conditions for the application of the solution (5). From 
its definition, @(t) can be written 





(a +i)(Q+18) rm _; 9-wO 
(ow +i)(Q,,+78,,) Q,+00, 5 Q 


m m 


(1 


m'| w,, 


on using (12). From this and equation (10) it is seen that the imaginary 
parts of ® on 7 = 0 and 7» = K’ are given respectively by 


im®, — 6*/(Q,,)y,  im@g. = 0; (16) 


while from (11) and (13) we have 
re(D,.~;-—®_,.) = X/Q),. (17) 


Equations (16) and (17) show that the solution (15) is applicable to ® once 
the function f 
To obtain f 


is determined. 


m 


consider the function 


In f,, 4 In(Q? +-62,)4 itan-a(o), 


m 


m? 


an analytic function of ¢. From (12) and (13) the boundary conditions for 
Inf , 


are: 


m 


on n , —2K y < 2K im Inf,, = 9, 
imInf,, = —tan~!a, 
imlnf,, = tan-a; 


re[ (Inf, )ox— (In fin) ox] = 0. 
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Hence a solution may be obtained from equation (5). We find 
0 


1-+-en(y*—t+ik’)dn(y*—t+-71K’) zs 
sn(y*—t+-1K’) 





inf,,(t) = A,,-4 tan-‘e@ [| 
2a : 
2K 
—n ] 1 
+Z(y*—t4+1K we... dy* —— tan- aw x 
2K 27 


2K 


* [l+en(y*—t+iK’)dn(y*—t+ik’) 
7 sn(y*—t+ik’) 

0 

After performing the necessary integrations we obtain 


l : 1—kedt 
i i. memeeaemenane Be 8 
Inf, = Ay + tan win) (18) 





yv 
A 


4 By*—t4iK') + a4". 


If we define a quantity « by 
tan-'w = 7e, (19) 
. 1—kedt\* 
then (18) gives Ahh a mgraen hg 
09) 8 Jn (Teer) 

where A, = expA,,. 

We are now in a position to calculate the values of f,, on the boundaries, 
and hence, from (16) and (17), a formula for ®(¢). On 

n= 0,-—-2K <y < 2K, 

we have @,, = 0, and so from (20), 


ay 1—kedy\* 
Om = A(isFea) 
1 * 


5 € 
Thus by (16), im®, = y My (az) ; 
Similarly on y = +2K, 0 < yn < K’, we have 


aad 1-+-ked in\* 
ee 4 poar) 


(20) 


so that by (17), 
1 1—kedin\* 
®@ -—QO = _—_ a—__—_—__" . 22 
re(Prx—O-2x) a*(pare) 
Using equations (21) and (22), together with im®,. = 0, we find from (5) 
that a solution for ® is 





(| +ked a) [! +en(y* —t)dn(y* —t) 


+2(y*—t)|dy*+ 


1—kedy* sn(y*—t) 





_(1—ked in*\£ [en(in* —t)dn(in* —t)—1 bye 
~ \l+kedin* sn(in*—t) 
_en(in* +t)dn(in* +t)—1 
sn(t7n* +-t) 


. 





+ Hn) Rint +0]. (23) 
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This equation, together with (14) and (20), immediately yields an ex- 
pression for f(t). After rearranging the terms with the aid of the addition 
properties of the elliptic functions we obtain 





2K 
f(t) ( ape 1 [ o(! Lked y*\* en y* dn y* Lentdnt 
i l+kedt} | 2m . l—ked y* sn y*—snt 

2K 


~ Z(t) +-Z(y*) |dy* 4 





K? 
1 f ,f/l—kedtin*\€ entdnt—enin* dnin* 
t= | x( l snt et ] 


Z(t)\d «| (24 
l+kedin*, (| 7 





sn*i7n* —sn*t 
0 


(where A, AA,,). 


m 


4. Conditions on the function /(/) 

We first apply the condition that the flow upstream is undisturbed. 
Taking the limit ¢ + 7K’ in equation (24), and equating real and in» zinary 
parts, we find 





2K 
l E a(i- kedy*\* sa . m i 
Ay+— | 0 | a) [—ksny*+Zy*)] dy* = 0 (28) 
Zn 1—kedy 
2K 
2K ; K’ , ‘ 
w | bad «,*\€ . - . : € 
and | a*(_ bed) dy*4 | x(; bed sy dy* 0. (26) 
} l ked y* ; 14 kedin* 
2K 0 
With the aid of these equations (24) can be written 
f()= 
n €| . | bad .,* 1 ay * * la r 
=f =a) | or: kedy ) eny dn y*. entent . seny* dy*4 
2r\l+kedt/| . l—ked y* sn y*—snft ‘ 
IK 
2 pm ; 
a baad font ; a ; 
| x(' kedtn sntcntdne ~en ty dnin dy*\. 
J 1+-kedin* sn*in* —sn*t 
0 
(27) 


Further conditions on f(t) may be obtained by considering the circula- 
tion. If the unsteady motion is assumed to start from rest at time ¢t = 0, 
then the circulation about the aerofoil is initially zero. At any subsequent 
time ¢ the circulation about any contour enclosing the aerofoil and (finite) 
wake will, by Kelvin’s theorem, also be zero. In particular, if the unsteady 
motion has persisted for an indefinitely long time, the contour vill extend 
downstream to infinity. Hence we have 

( dw = 0, 
c 
where w is the complex potential function, and C is the closed circuit 
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containing the aerofoil and vortex sheet, shown in Fig. 3(a). Since 
CG) O23 
—_ = _e~i8 4 eS, 
dz U 
therefore | etdz= 0. 
c 














(a) z plane 











Cc 


a 
<= 








a 





(b) tplane 
Fic. 3 
To first order in the perturbations the last equation may be written 


[ fe) dz = 0, 
Cc 


since { dz = 0. Differentiation of equation (6) gives 


Cc 
7 
= ——sntdt 
dz ya dt, (28) 
so that we have fre snt dt = 0, ‘ 
é 


where C is now the path shown in Fig. 3(b). Equating real and imaginary 
parts of this equation gives 


2K . 
| *sny dy = 0 (29) 
—2K 
2K K’ 
and | QAy)sny dy—i {| Xsnin dy = 0. (30) 
-2K 0 


5092 .46 M 
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The form of equation (30) is unsatisfactory as it contains the unknown 
function Q(y). Hence we substitute the appropriate limits from (27) into 
this — and express it as 


2K 


2h 


[6 ar: +-ked y*\€ [Gr en x enydny+en y* dn y* 
~ Feds?) 3) e a sn y*—sny 
2K 





dydy* 





2K 
[l—ked ony ( 1—k eat enty BI 11” dydn* 


1 +kedin* ltkedy sn*i7* —sn*y 
2K 


=z @. (31) 
This equation represents a further integral condition relating 0* and X. 
Finally, the value of X is subject to the condition that the pressure is 


continuous across the vortex sheet. From Bernoulli’s equation we have 
that on either side of the sheet 


l cp 


0 (y= 0,2 > 2a), 
pcx 


so that across the sheet we obtain, to first order, 


eX 1 eX | 


2 0. 
ex v2 et 


In terms of the t-plane variables, this is, by equation (28), 


ins? oX , ox 0 
s —% = V. 
The 1 on ot 


5. Pressure on the aerofoil surface 


For an upwash of general form, an exact expression may now be derived 
for the pressure distribution on the aerofoil surface. Bernoulli’s equation 
can be written 


— of 
p = 4p(U*—q")—p= | ads, (34) 
the constant —* omitted as it makes no contribution to the lift and 
moment. Writing s = x on the aerofoil surface, we obtain from equation 


(28) 


ds = dx = = sn y dy, 
7 


‘ 


and so to first order in the perturbation i q—U, (34) can be written 


p = pURFQ+ = any, (35) 
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2K 
- 1+ked y*\*/1—kedy ‘jen ydn y+en y* dn y* | 
fel. Mkt! r 4 ° 
6 ( cag (i= per) | sny*—sny ksny ju" + 


K? 
1 [2 (Spears) (eet) en y dn y—enin* dnin* dn*. 


T Se | 1+kedin* sn*i7n* —sn*y 
0 


1+kedy 





(36) 


We now substitute for Q from (36) into (35). It is shown in Appendix 2 
that an expression for the pressure distribution is 


9K 


0 Tf .,,(1+kedy*\*/l1—kedy\€ 
p 1 Y at 
» < - l 6* +G* 

F 2 | \ (i peal (pe) " 


=7T 





9K 





F jn y+en y* dn y* 


a bsn y*|dy*+ 
sn y*—sny war poe \y 





: 3K x 
_ pUck [ (U0*+G*) pera) | (eae = 
7 1—kedy* l+kedy 


x fem ydny-+sn 72 y)—- Prey" (37) 


y* 
» le 
where we have put G(y*) = ahh [ 6* sn y dy, (38) 
7 
0 
and where F, H are constants defined by 


* /1—ked 

—kC y € 
| (en) = 
i 


and =) [enydny+sn yZ(y)| dy. (40) 
2K 


It is an important feature that equation (37) for p contains no terms 
involving X, and so the expressions for the lift and moment will similarly 
be free of such terms. 


6. Movement of the front stagnation point 

It is clear from equation (37) that the pressure, and hence the lift and 
moment, are determined once the flow direction at the aerofoil surface, 
6*, is given. However 6* is not completely known a priori, since it is 
composed not only of a term arising from the prescribed upwash velocity, 
but also of a term associated with the motion of the front stagnation point. 
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The latter contribution is not prescribed, and must be calculated from 
the conditions of the flow. 

The position of the front stagnation point in the mean steady flow is 
at y* = 0 in the ¢ plane. Suppose that at any instant of the unsteady 
motion its position is at y* = —6, where 6 will be small, since the ampli- 
tude of the unsteady motion is small. Then in the range —5 < y* < 0, 
the flow direction will be reversed, that is, 6* will be increased by z in 
this range. 

if v(y*,t) be the upwash distribution, then apart from the above term 
the flow direction at the aerofoil is 6* = v/U to first order. Taking a 
harmonic upwash of the form 


v(y, t) vo(y)e™, 
then at the aerofoil surface we may write 
(Yeo (2K <y <2K) 
O* —= l fe 
& (—8 <y < 0). 


Moreover equations (32) and (33) enable us to take X in the form 


X= Xo int—2,0) — X,e( “" ’ (42) 
1+-kedin 
where X, is a constant. 
The quantity 4 can be uniquely determined from equations (26) and 
(31). It is found convenient to replace (31) by a linear combination of it 
with (26), namely, 


2K 


{ g+(Ltkedy*\< 
(i —ked s 


- i 
~2K —2 





(; ba) | 2 enydny+eny*dny* | 
8 + 


l+kedy) | sn y*—sny 
+enydn 7| dydy* +- 
K’ 2K 
f y(1—kedin*\* f (1—kedy\< ” 
* | 1-+kedin* | lt+kedy) ~ 
0 —2K 


aan Se * a 
- enty BY enin* dnin 





sn%i7*—an?y +enydn 7| dydy* = 0. (43) 
To evaluate the integrals occurring here, as well as the constants F and 
H, it is necessary to expand the integrands in series and integrate term 
by term. Since thickness effects are being neglected, we may assume our 
results to be valid for small values of the ratio aerofoil length to tunnel 
width. Hence we expand in powers of the elliptic parameter k, defined by 
equation (7) to be a direct function of this ratio. The basic approximations 
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are found to be 


L—kedu\* _ | _2¢ked u+2%k* ed?u+O(i?) ates 
l+kedu 


and (oe ' - tL %ekedu+2 ek? ed?u+ O(k*). (45) 
l—kedu 


The formulae (41) and (42) for 6* and X are now substituted into 
equations (26) and (41). Using (44) and (45) we carry out the integrations 
and obtain, to first order in 8, 

a Ok 
+ 22h) TH | vo(y*)[ 1+ 2ek en y* — 2e?k? sn*y* | dy*+- 


. 


2K 


K’ 
7 l1—kedi €+Hivh/rU) 
1 ¥ iv 7 — 0 

Xve | (spar) ” (46) 


2K 
pivt e 


2e*k?) 7 | Vo(y*)[en y* dn y*+-sn y* Z(y*) + 
2K 
+- 2ek( 1 — 2e*k* sn*y*)] dy* — 


K’ 
1 . , fe =" di e+ivh/nl’) 
F Xe" | = 2 x 


1+kedin 





l—kedy >» enydny—enin* dnin* 
dn y| dydyn* = 
(= " rope) [sn*y sn*in* —sn*y acted Vad 
(47) 
By elimination of X,, 5 may be obtained from these equations. 
The integrals involving 7 may be evaluated exactly in terms of Legendre 
functions. It is shown in Appendix 3 that 


K’ 2K 

[ (= 7 kod in\eHonn0) { eat) x 
J} \l+kedin l+kedy, 
0 ~2K 


enydny—ecnin* dnin*® , 
x fanty ; antiq*—anty ‘ +enydn7} wr 





[enindnin+sninZ(in)| dn— 


1—ked in\eH™” 
-5 a | (see a, 


where we have put r. (48) 


1—kedin\tH") 
1+kedin 
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Transforming to the variable £, defined by 
x 

aa * 

we find K’ ? 
1—kedin €+(id/r) r e-fer+id dé 

na | (ibedin) gy of eterna 
1+kedin v2k’ J ,/{cosh r§—cosh r} 
i 
and 


os an €+iA/r 
= aie [enindnin-+sninZ(in)|dy 
: dir 


0 
¢ &er ayt sinh r€ 


rk’ v20 if’... 7E\]) 
: ' be nccrnend 4+.—— Z| en-{ — sinh = : 
2V2k f \,{eosh r€é—coshr}  k’ zen (7 a *)| 7 
i 


It has been shown (5) that these integrals are expressible as 


i, = 7 Qx(cosh r) 


and 
rk’ , 2E , 
I, ia =l Qy_,(cosh r)+4 (ga 1} Qs (cosh r)—(N + 1)Q,(cosh nl. 


where (,(z) is the Legendre function of the second kind, Z is the complete 
elliptic integral of the second kind, and 
iA 


N te—4. (49) 
r 
If for brevity we deffhe eK 
; ~ 
A i | Vol 1 —2€2k? sn*y] dy (50) 
E 2K 
2K 
] ° 
and b, 7 | volen y dn y+sn yZ(y)] dy, (51) 
7 PA 
then equations (46) and (47) may now be written, correct to order k?, 
o(1+k\* Xe 1 
a( —— +-(1-+-2e*k*)e™ (a, + 2ekb,) + Kee”. Oy 0 (52) 
1—k; an ok 


and 


a(; x i) +(1-+ 2e2k* be Midy T 2ekay) = 


eee tivo, ton cH > | 
= nae Qyat (aan? Jex— (N+1)Qy. | — TF Qn} = 0. 


(53) 
Eliminating X, from (52) and (53) yields, to the required order of accuracy, 
$= — e™ {4 (a,+ by)+ 4(ay—b,)T(A, r,e)}, (54) 
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where 





‘ ~ 2E 4ekH 4idrk 
Tir a ‘) 1—2ek N Qya—(A ‘a NQvat(eea— 1 PE? =F) @s 
sis 1+2k — . 2E 4ckH  4iAk\,._ 
NQyi—(N+1 NOvat (Re —1— Fr? + ra) Qy 
(55) 


7. The lift and moment 
The lift Z, and moment M, about the mid-chord point are given by 
= x 
L=— | p dx(y), M = | xp dx(y). 
y=-2K y= —2K 
Using equation (28) the lift may be written 


ba 


7 


ad 
| psny dy. (56) 
2K 


By using the expansions (44) and (45) in equation (37) we obtain an 
approximate expression for the pressure p. Substituting this expression 
into (56) we get eventually 








2K 
L= Saal | (U0*+G*)[en y* dn y*-+sn y*Z(y*)—ek? ed y*] dy*, 
wn 
—2K 
(57) 
which is the result correct to order k*. Similarly the moment can, by 
(28), be written 2K 
(h\? ¢ 1—ked 
m=2(*) | vim 
;) | psn yin( 79") dy (58) 
—2K 
which gives, after substituting for p, 
Khk%pU f 
4Kh*k*pU ’ , P 
M = — = | (Ud*+G*)[ 1 —2sn*y*+k*(}—} sn*y*)] dy*. 
—2K 
(59) 


For harmonic oscillations we may substitute for 6* and G* from (41). 
We find 


ae 4Kh 





bella ety 


7 
2K 


(vo +go)[en y* dn y* +-sn y*Z(y*)—ek? ed y*)ay}, (60) 
—~2K 


e ivt 


+o. 
al 
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and 
*p21.2, 772 
M = ar 51-4 4k2) 4 
4 
2K 
Ff La 9 ap 2ay* 1 Ee2 4,*) *| ; 
+ | Wot dell1—2sn%y*+kd—dsnty*)] dy*!, (61) 
~2K 
y* 
Oo. 25 P 
where now Go as | vysn y dy. (62) 
nU J 
0 


8. Rigid-body oscillations 

As a particular case we consider a rigid aerofoil undergoing a harmonic 
oscillation compounded of a flapping motion parallel to the z-axis and 
a pitching about the mid-chord point. That is, 


y p ye i vt x = we ivt 


respectively, where ¥°, «° are the amplitudes. The upwash velocity in this 
case is obtained as follows: The velocity normal to the surface of a point 
P on the surface distant x from the mid-point is, to first order, 


Yy—xa. 


Again, the fluid velocity normal to the surface at this point is, to the 


same order, Vex 
X-+ U. 


Equating these expressions, 
v — [7 x%ivt 4 iv(y®—xa°)e™, 
That is, 


Vp = tro(y®—zra°)—Uo® = in| yo“ In a er)|— YU. (63) 


Substituting from (63) into (50), (51), and (62) gives 


4K __ two hk = 


ay = (ivy°—o°V) —— (1—ek*), sb 








al 
; Ws 5 oss 1—kedy\ , v*h?a® 1—ked y\]}? 
= — (ivy’—o°l ot A Ntlcclaall eantaoel 4 | 
We also have the expansion (8) 
= $n[1+}h2+ 0(k4)] 
and h =“ [1— 42+ O(K)] 


k 


which is derived from equation (7). Using these results in (60) and (61) 
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we eventually obtain 


L = —napUe'™| 2(ivy—a®U)(1+-T)(1-+ 4k? —2e2k*) + 
4 2ivaa1—T)(1+ Ak? ek?) + te (ivy®—20°U)(1-+,k2—e2k2)| (64) 
and 
M = —na%U* i (ivy? —a°U)(1+- T)(1-+ yk? —e2k2) + 
2ivan 2v?a?a°| e 
odie, conaguasame _ pee 6 
(1—T) vu | (65) 


We express the lift and moment more conveniently as non-dimensional 
air-load coefticients. The flapping motion is described by the coefficients 
L,, lo, my, m, and the pitching by Js, 1,, “ m,, defined by 

L 
2napU [jigivt ~ 

M 
ina®pU%m 


= (1, +0, Ey +il,)a? , 


= (on, img) E+ (rmg+img)a?. (66) 
Replacing v by the dimensionless frequency parameter A = 2av/U, we find 
from (64)—(66), 

1, = Aim T+22+ Ak?A(2im T'+A)(1— 1 2¢2) 
1, = —A(1-+-re 7) —3k2A(1-+re 7')(1—12e2) 
l, = 1+re a T+ Ak?(1—12c*)(24+2 re T—4Aim 7) 
1, = im T+ 3A(3+1e T)+ 4 k?(1— 12e*)[ 2 im T +-4A(3+re T)| 
m, = sAim T+-3k?(1—12e*)Aim T 
m, = —4A(1+re the jjk*(1— 12€?)A(1+-re T’) 
ms = 4(1+reT)—}Aim T+ JA?2+ A k*(1— 12€?)(1+ re 7) 
m, = tim T— hig T)+-3,k7(1—12¢*)im T } 


On putting « = 0, the equations (67) immediately reduce to the corre- 
sponding formulae for a solid wall tunnel (5). 


. (67) 





APPENDIX 1 


TRANSFORMATION OF THE SOLUTION Fit) 
A function which satisfies, in the rectangle —2K < y < 2K, 0 < » < K’ the 
conditions described in section 3 above has been found (7) to be, in our notation, 


F(t) A+ap £ (im Fe x” —t)\n]- im Fu [7k = (y* —1)|r]} ay + 


a re(Fex—F ax) (3![ Ze (in* —t)7,]-# 5, Kin" —t)|r]} dye , (68) 


0 
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where the theta functions are as defined in (9), p. 463, and have parameter 
K 

2K" 


T 


This solution can be transformed into an alternative form by using the properties 
of the theta functions given in (9), chap. 21. From Landen’s transformation we 
have 


3, (ul7,).3,(ul7,) = const. x F,(2u|r), 
rr , . e..- ° 
where 7 = 27, i‘K’ /K. Differentiating logarithmically, 
wy H w wv, 
ult, )+_q (ult) 2 5- (2u |r) 21+ 2 = (2u+ §nr}r) 
5, ' ae * : 3, 3, . 
21 + 293(0!r)Z[93(0\7)(2u + }r7)], (69) 


from the definition of the zeta function. An equation involving theta functions (9) 
(loe. cit.) is 


Q’ Qc ar 
o 9 3 F,(y (2) (y + 2) 
=*(y) + = (z)— = (y+z) = 8,(0)9,(0) 2 sd i< 
vy; Vv; Vv, F,(y) a (zy, (y 
If we put y = u+ rr, z kv and use the properties of the theta functions, this 
becomes 
o x H.(ulr, 9 (uitr,) \3 {O 
l 2 a2 3 ‘1 4 ‘1 2 3\ r Qo 
a (téiT;) (uiT,) v3(OirT,) —— P(Olr 5(930 T)2u}, 
a, : o, ; : : H(t) 7, g(t |7,) » 1) 3. (O\7r . 


(70) 


from Landen’s transformation and the definition of the elliptic function. It is found 


that fe 2K oa %(O\r) 4K 
v;3(0\7) . vi(0\7,) ; 
7 . 3,(0\7) 7 


further, the periodicity of the zeta function gives 
: . , SRE en 2udn 2u itr 
Z(2u+7kK) Z(2u+ikK’) Z(2u)4 a 


sn2u 2K - 








Substituting these results into (69) and (70), we obtain 











4Ku 4Ku 
cy Ku) cn -— dn —— +1 
—— (u\T,) Z\— = (71) 
a, a 4Ku | 
sn - 
4Ku 4Ku , 
- ey 2K (tk) cn = — | = 
an u't,) / ° (72 
i 3, 1 | Ku | ) 
sn —— 
By a similar procedure it can be shown that 
ae One. wa 
5 oK iKu en| = iK }an( - tik ) +] al 
a (u|t;) Z| LiK’)+ . ik yak (73) 
i 7 7 a 2 
' sn| 4 ik’) 


7 


The theta functions of equation (68) may now be replaced by the expressions (71) 
(73), and this leads immediately to the form (5). 
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APPENDIX 2 
AN EXPRESSION FOR THE PRESSURE 


For the sake of brevity we define 














1—kedu)\‘(1—ked v\f(snu(envdnv—enudnu) 
filu,v) | 1+ tain) (; + =) \ sn@u—sn*v i 
1+kedu\*(1-—-kedv\‘(envdne+cnudnu : | 
falu, v) 7 teas) (; a) shu—snv +h +Sene; 
and 
l+kedu\‘/1—kedv\*{envdnv+enudnu | 
f.(u,v) (j Seda) ( + Fda) \ snu—snv — Biv) + ben ~ 
Then with some rearrangement rr (36) may be written 
2K 
l r sa 1+kedy* —kedy\* 
) * ,* r, * i‘ 
" 2a O°faly*.y) dy® | 6*( “rei (r kedy ] ey 


kK’ 





x(u ~kedin* ) (ss) 


l l 
4 F ieee healt ee 
2a X sn in*nsyfilin®.y)dq®4 2a g 1+kedin* kedy 


0 
< fenin* dnin* +snin*Z(in*)—enydny]dn*. (74) 


From equation (26) it is clear that 





enydny [ x(': l ~ked in* \'( ae @ 





2arsny , + ked in* kedy 
2K 
= enydny ofp teer (EG ESst) ay o 
2rsny 1—kedy* kedy . 
while rearrangement of (31) enables us to write 
x 
] 5 x l kedin* ~-kedy se —_ * 
sar) (5 Seas) ( tkedy J" om ing® dn in +0 9249) dy 
0 
hose fF ts eee . 
“SFr kedy’ , mediate yee Ot 
0 og 
l 1—kedy' ./ 1+kedy* J . 
Se rhea) | 0 \ B( Fea ) + anu faly*,u)duid 
Using these results in (74) we find i 
a raf 1—kedy\« f \ 
- Cae = *e * ~7o) | * 1 





1+kedy* Gen a) i 
I ‘e(! kedy* :) I+kedy enydny+snyZ(y)—, dy a 


ee i sn in*f(in®, u)dul dy? 
Kk 


(75) 


— * 
ct ryeninthalin®, n— Wipe) 


0 -2 
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Now equation (33) gives that 
eX 
ct 


nU 
~ Qhk 


insin 


and hence the last integral of (75) can be written 


me 
I Bh [ OX | * 4) . <==! s) 
1 SRR JS ay8 TOY PATS kody) 
0 
~ 
j { Fi. kedy 
; F\i 


€ 
X —= fi(in*.y)— 
4hk | | on* it Bd tkedy 
0 


the integrated term vanishing. 
-kedin 


] 
pha “i 


<fenvdnv 


‘ 
—= fi(v, im) 2ek( 
cv 


Thus we obtain 


2K 


k 
I C | x2 fy, in*)an* fe 
‘ 7 y.4 YY + 

t thk . ey" i\Y¥ 7) 7 x 


. [enyany +snyZ(y) 


where the 


eancelled out. 
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, 


It can be verified by differentiation that 


\iFeay* 


cy 


2K 
| Silin*, u) du dy* 


4 
a 
§——- in*,u)du} dn*, 
2K 


kedv\* 


kedv 


sn vZ(v) sninZ(in)). 


en in dn in 


kedy)\‘ 


tkedy 


-kedy* 


) ( 


= dy*, (76) 


last term is found by replacement from (26), and other terms have 


From the above definitions it is readily seen that 





; tkedu)\* tend 
sn vf,(u, v) sn uf.(u, v) (; ee) '(: (! 7) 
<fenudnu+snuZ(u)+envdne+snvZ(v)}. 
This allows the first two terms of (75) to be written 
2K aif 2K 
kedy 
I, fo *\sn-y*faly (y*,y)— at =r) | rer ey 
2K 2K 
Introduce now the function 
y* 
2hk 
G(y*) — O* sn y* dy*; 
“ 
then we can put 
I f Slew J ens | by*, \dul dy 
. —\ f(y*,y)——\- uu) ¢ 
” ii x dy* |: Peel FN kedy fly —— x 
1 ye! _é 11 kody): T . \ 
- Thk F lay = Sly*.y)— Ar red y) dé ss Saly*, pene ays 
~2K 


the integrated term vanishing. 
¢ l+ked 
—%, Saly*.»)— 2er(— oer * 


—kedy* 


bal 


fily*,v) = 


=/2 


x [eny* dn y* 





Again it can be 


) (5 


shown that 


—k sy" 


+kedv 


+sny*Z(y*)+envdnv+snvZ(v)]. 
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Hence we find 
2K 





2 


| ‘ { 1+kedy*\</1— ant) 
em ye) * 1 OL 
thk | Og Sly) 2ek (Eo) \>Feay 
2K 


x [enyany +snyZ(y)—“| dy*. (77) 


From (76) and (77) we obtain 














y x’ 2K 
2QhkpU [ Q U 
"se any eo | Xfily,in*) dy e | G*fsly*,y)dy* + 
0 é te 
2K y 
ekpU [ l+kedy*)\* [ haat a) 
a na T * a _H °. 
eg w | ( 0 +G i teas) | (= kedy enydny+snyZ(y PF dydy 
—2K 


Adding this to pU*Q, given from (36), yields the required formula for the pressure. 


APPENDIX 3 
EVALUATION OF AN INTEGRAL 


Let 


2K 


K’ OK 
F —kedin\+* f (1 —— enydny—enindn in 
ad (+ “Eediy) | ¢ tkedy sn*y sn? in —sn*y +enydny/dydn, 





0 
where a = ivh/7U. This may be written 


a 
Kus 1+H | ((—s3) dy, 


1+ kedin 
0 
where K’ on 
* (l—kedin\tt= [ fe 
’ | (; +ked i) | (= kedy sn y gly, in) dydn 
. —2K 
and gis, 0) matt vdnv—cnudnu ~ Zu). 





sn*u—sn*v 
t (i-eten ‘rr Ate} 
ink iFbae) _ mms 1+kedy vty 





Now 





K’ 2K 
cB 1 { f[ (l—kedin “(37)'< 
~ ~~ Dhe i | (; Eoin) rp Bade) 5 By My in) dydn, 


0 —2K 
the integrated term vanishing. We can show by differentiation that 


< 5 MY in) = iz inn 





KR’ 2K 
ue ee , (1—kedin “(= =2y'= 
mice dis: ax | | (Fare) 1+kedy on my) dydn 


0 —2Kk 


( F (1—kedins+*/1—ked y) 
=<] | (ey (a2) smn atin 9) dyn. 
0 





€ 1+kedin 1+kedy 





174 THE AERODYNAMIC FORCES ON AN AEROFOIL 


From the definition of g we have 
sn in g(ty, y) sny g(y, in) +enydny+snyZ(y)—en in dn in—sninZ(in), 
so that 


= 

€+a@ r fe kedin €+o 
I= — E -H | (— bedie) dy 
0 


K’ 
° (1—kedin\**" : . 
F | { — Todi) {en in dnin+sninZ(in)} dy |. 
0 
This gives 
K K’ 
€ _{ (l—kedin\*:, ; : ene ; fa 
F | (3 Fedin) [enin dnin+sninZ(in)|dn—H | ( 


0 0 


ked 1) . 
t ked ‘7 


an}. 


which leads directly to the required expression for J. 
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SUMMARY 

Exact solutions of the Navier-Stokes equations are derived for two-dimensional 
flow against an infinite flat plate normal to the stream and performing an arbitrary 
transverse motion. This generalizes Glauert’s solutioa for an oscillating transvers- 
motion. 

After a brief approximate investigation of the reaction of an arbitrary boundary 
layer to arbitrary changes in the free-stream velocity, the problem described above 
is considered. When the wall moves from rest, then, by means of a Laplace-trans- 
form technique, expansions of the velocity distribution for small and large times 
are given in terms of the velocity of the wall. Further, it is indicated how a Pohl- 
hausen type of method may always be used to obtain an approximate solution for 
the purpose of linking up these expansions across the range of times for which neither 
is valid. The method may also be used for motions other than those starting from 
rest. In uniform motion started impulsively from rest, the expansions overlap and 
the approximate solution agrees well with them. 


1. Introduction 

For many years the study of unsteady laminar bovndary layers was 
restricted to problems of boundary-layer growth in motion from rest and 
the boundary layer in oscillatory motions with zero mean velocity. 
tecently Lighthill (1), Stuart (2), Glauert (3), and Rott (4) have investi- 
gated how the boundary layer on a two-dimensional body reacts to 
fluctuations in the main stream. Glauert (3) and Rott (4) in particular 
studied the two-dimensional flow when an infinite plane wall normal to 
the main stream makes transverse oscillations in its own plane. 

Watson (5) has generalized the solution given by Stuart (2) to the case 
of a main-stream flow which varies arbitrarily with time. In the case of 
an impulsive increase of the main-stream velocity from one value to 
another, it was found that the boundary layer initially reacts as though 
it were inviscid, except within a narrow region adjacent to the surface. 
Here a secondary boundary layer of the classic Rayleigh impulsive type 
is formed. As time progresses the secondary layer increases in thickness 
and gradually distorts the flow until the ultimate steady-state boundary 
layer is attained. 

The primary purpose of this paper is to generalize the work of Glauert 

[Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 2, 1959) 
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(3) and Rott (4) to the case of arbitrary motions of the wall. However, 
a preliminary general analysis of the unsteady two-dimensional laminar 
boundary-layer equations is given in section 2. For an arbitrary motion 
of a body commencing at t = 0, the first term in a series expansion for 
the velocity in the boundary layer at small times is obtained by omitting 
the convection terms. Further, for the particular case of a small impulsive 
increase in free-stream velocity the first term in an expansion for large 
times is derived from the quasi-steady solution given by Lighthill. 

The generalization of Glauert’s problem is considered in section 3, where 
a linear partial differential equation is obtained. In section 4, by using 
the Laplace-transform technique, the velocity, as a series expansion for 
small times and as an asymptotic expansion for large times, is easily 
derived from Glauert’s solutions for oscillations of small and large fre- 
quencies. Numerical results are given for the case of uniform motion of 
the wall, started impulsively from rest. Examination of the curve of skin 
friction against time shows that the expansion for small times tends 
smoothly to the constant limiting value for large times. 

In general, however, we cannot expect the two solutions to join quite 
so conveniently, for which reason an alternative approximate solution valid 
at all times is given in section 5. This is a Pohlhausen type of method 
in which the differential equation describing the unsteady flow is satisfied 
only at the wall, at infinity, and on an average by integrating it across 
the boundary layer. This approximate solution, being valid at all times, 
may be used to interpolate accurately between the two expansions when- 
ever there is an interval of time for which neither is applicable. 

For impulsive motion of the wall the curves of skin friction given by 
the exact (Laplace) and approximate (Pohlhausen) methods are near to 
each other and, to a good approximation, are parallel. This encourages 
us to hope that in other cases reasonably easy interpolation will be possible, 
particularly as a similar result was obtained by Glauert and Lighthill (6) 
in using a Pohlhausen method to join two expansions for the axisym- 
metric boundary layer on a long circular cylinder. For impulsive motion 
of the wall the characteristic time t, representative of the time taken for 
the boundary layer to settle down to its steady limiting state is of order 
53,/v, where 5,, is the undisturbed displacement thickness. A similar result 
was found for the asymptotic suction profile (Watson (5)). 


2. An analysis of the unsteady two-dimensional laminar boundary 
layer equations 
Let x denote the distance along a body, z the distance normal to it, 
u and w the corresponding components of velocity, t the time, p the 
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pressure, p the density, and v the kinematic viscosity. For incompressible 
flow, the equations of the boundary layer and of continuity are 

pm U7, aw 


me wo ae 
oe et Oe ea 


' 2.1 
a (2.1) 


(2.2) 
subject to the boundary conditions 


u=-w=—0 at z=0, u—> U(z,t) as z>o, (2.3) 
where LU (x,t) denotes the velocity at the edge of the boundary layer. 
We now write 
u(x,2,t) = Uo(x,z)+-u, (2, z, t) 
w(x,z,t) = we(x,z)+w,(2,2z,t) }> (2.4) 
U(x,t) = U(x) +U,(2, t) 
where the suffix 0 denotes a basic or original steady motion, and the suffix 


| denotes the time-dependent part of the motion. The steady part of the 
motion satisfies the equations 


dUy | up 


"Fe 


subject to 


Ho=w=0 at 3 = 0, Uy > U(x) as 22>. (2.6) 


When equations (2.4) are substituted into (2.1) to (2.3), and (2.5) and (2.6) 
are used, we find that 


( ou ou, dU, , ou, 
t u,—°+ w,—! + w, - e__4_Uj—t— 
Oz cz at Ox 
as U U, . mes 
= {32> =< 
Ox ox 


ou ou 
ee WT oe 2 


ot ° Ox ox 


_- = 0, 
dx 2 
subject to the boundary conditions 
u=w=0 atz=0, u, > U(z,t) asz>o. (2.9) 
We now suppose that the flow is steady for t < 0, so that the time- 
dependent motion commences at ¢ = 0. Additional vorticity will be 
created at the surface of the body and will diffuse outwards, reacting with 
the main boundary layer. This region of additional vorticity created by 


5092 .46 N 
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the time-dependent part of the motion will be termed the ‘secondary 
boundary layer’, outside of which the flow reacts to the time-dependent 
part of the motion as if inviscid; accordingly, in this region u, = U,(z, t). 
For ¢ small the secondary layer will be very thin and the viscous term in 
(2.7) very large compared with the convection terms. Thus for such values 
of ¢ equation (2.7) may be replaced by 
éu, eu, eu, 


== Jb 


z>0,t > 0), 
ot ot 0z* \ ) 


which may be written in the form 


subject to U,(z,t) when 


as > 0, ’> ¢, 
when z > 0,t = 0. 


The solution of (2.10) (ef. (7) p. 334) which satisfies these boundary condi- 
tions is t 


. 


ul, U(7,)—s | (2, t— Ir. (2.11) 


0 
From (2.8), (2.9), and (2.11) it can be shown that 
f 


Cx 


T. 
7 


Cx 7 
0 


Provided that ¢ is small enough for equation (2.10) to be valid, the total 
velocity components are given by (2.4), (2.11), and (2.12). 


_ 


It follows from (2.4) that the skin friction is given by 


_ (=) (2.13) 
CZ 2-0 


where 7,,, is the steady skin-friction. To evaluate (2.13) for the particular 
solution (2.11), it is sufficient to assume the mild condition 


U(x, t)- U(x, t -r)|} < K\ir|*, 
where K, k are positive constants and k > 3; then it can readily be shown 
that the skin friction, for small and positive times, is given by 


t 


Cen, i f 


(mt)! 2(av)! 
0 


(U,(2, 1) —U x, —2)} (2.14) 


(T~— Ti) p= 
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Consider now the particular case in which the free-stream velocity in- 
creases impulsively by AU,(x) at t = 0 and remains constant so that 


U,(x,t) = 0 for t < 0, 
U(x, t) = AU (x) fort > 0. 

In this case (2.11) yields 
u, = AU (x)erf{z/2(vt)t}, 


a 


where erfx = 22-4 | e- dx. 
0 

Similarly, by using (2.12) or by using (2.8) directly, we can obtain a 
formula for the velocity component w,. The skin friction is determined 
immediately from (2.14) as (7,,—7,,,)/u = AU,(x)(zvt)-* for small times. 
This secondary boundary layer is, of course, the Rayleigh impulsive flow; 
the theory shows that, initially, it may be added linearly to the original 
steady flow. 

For large values of the time the velocity distribution will tend to that 
for steady flow with free-stream velocity 

U(ax,t) = (14+-A)U,(2). 

If \ is small compared with unity then u,, w, will be small compared 
with up, wy. On neglecting the squares and products of u,, w,, U, in (2.7) 
we deduce that 

Uy (0 em), w, = 48(wo+22%), 

oz oz 

the ‘quasi-steady’ solution obtained by Lighthill (1). The skin friction is 
found to be 7,, = 7,,,+3Ar,,,/2. 

The complete solution of equations (2.7) and (2.8) for all times presents 
considerable difficulties even in the case of impulsive flow. A simple 
particular case, that of flow near a stagnation point when the wall is 
impulsively set into motion with a steady velocity, is solved in the follow- 
ing sections; this is equivalent to the case of U,(x) a x with U, = 0 for 
‘<0 and U, = constant for t > 0. The method is applicable to other 
cases of unsteady motion of the wall corresponding to the same U,(z), 
but with U\(2,t) as a function of time only. Because equation (2.7) is 
linear in these cases, exact solutions can be obtained; furthermore the 
solutions also satisfy the Navier-Stokes equations exactly. 


3. A general analysis of the flow near a stagnation point when the 
wall is in unsteady motion in its own plane 
The solution for the steady boundary layer near a stagnation point on 
an infinite plate normal to the stream has long been known, and recently 





180 J. WATSON 


Glauert (3) and Rott (4) have considered the unsteady flow which arises 
when the wall oscillates in its own plane. The solution, like that for the 
steady flow, is an exact solution of the Navier-Stokes equations. We will 
consider here a generalization of this oscillatory solution to the case in 
which the wall moves in an arbitrary manner. 


With coordinates and velocity components defined as in section 2, a 
steady solution of the Navier-Stokes and continuity equations is known 
in the form 

U, = az, Uy = axd'(n), We ~(av)*d(n), z(a/v)', (3.1) 
where a is a constant and ¢ satisfies 

d” 1 bd” | l —d’? 

subject to the conditions ; (3.2) 
$(0) = $'(0) = 0; (x) 

The corresponding pressure distribution (p) can also be calculated and is 


given by P—Po ~4p[a2x?+-av(?-+ 24’), (3.3) 


where p is the density and p, the stagnation pressure. 


For the case of flow in which the wall has a velocity which is an arbitrary 
function of time, it is possible to find an exact solution of the Navier— 
Stokes and continuity uations relative to axes fixed in space in the form 


u = axrd'(n)+eF(n,7), w —(av)*d(n), r= al, (3.4) 
where ¢ is a reference velocity for the wail and F is given by the equation 


F’1+46F’—¢'F ee 0. (3.5) 
CT 


Primes denote differentiation with respect to 7 and the boundary condi- 
Wonsare FF = Fils) at» = 0; F+0 asy>o, (3.6) 
where «F,(r) is the velocity of the wall. The pressure distribution is still 
given by equation (3.3). When (3.5), (3.6) are solved in conjunction with 
(3.2), the solution for arbitrary motion of the wall is obtained. 

Glauert (3) has investigated the case when the wall oscillates har- 
monically in which case (3.5) is separable, leaving an ordinary differential 
equation to be integrated with » as the independent variable. This same 
differential equation of course arises when the Fourier or Laplace transform 
is applied to (3.5). The latter will be used in this paper. 
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4. Solution by Laplace-transform technique 


We suppose that the motion of the wall starts at ¢ (or r) = 0. The 
Laplace transform, #(p), of the function x(r) is defined (7) to be 


ip) = | e?*2(7) dr,  re(p) >0. (4.1) 
0 


The inverse of this is 
yt+ia 


xr) = Orz(A) dX (y >e > 0). (4.2) 
mr 
y-iw 
On multiplying (3.5) by e-?7 and integrating with respect to 7 from 0 to 
oo we obtain F’+4F’—¢'F—pF = 0, (4.3) 
assuming that e-”’F(»,7) > 0 as r+ 00. From (3.6) the boundary condi- 
tions on F(», p) are 
F=FA(p) atn=0; F+0 asyn>o. (4.4) 
Accordingly the function ie ” 
F(n.p) = F(n, p)/Folp) (4.5) 
satisfies f’+¢4f'—¢f—pf = 0, (4.6) 
subject to the boundary conditions 
f=1 atyn=0 
f- 0 asyn>o 
A general solution of (4.6), (4.7) has not been obtained. However, an 
expansion of f for small and large p yields a formal expansion of F for 


small and large r. The solution for small 7 will be obtained in section 4.1, 
and that for large + in section 4.2. 


(4.7) 


4.1. Solution for small times 

If f(y, p) is expanded in a series for large p, and the resulting series is 
inverted term by term, then a formal solution of (3.5), (3.6) will be 
obtained for small positive t (Carslaw and Jaeger (7)). Now Glauert (3) 
has considered (4.6), (4.7) with p replaced by (iw/c). From his results we 


deduce that a solution of (4.6), (4.7) for large values of p can be obtained 
in the form 


2 n 
f(n,p) = 2. pie nat * Onn PY™y™ 
which can be rewritten 
F(n, Pp) = (> mm nm)e-ring (&, ammes 1”)p-te-Pin + 
i tJ 


+ > ( 2, Ammense 1)p-1-ine-vin (4.8) 
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with a change in the definition of n. The coefficients a,,,, are given in 
" 7 . 4 ; 4 
Table 1. We note that the inverses of e-?°", p-te-P*), p-l-ine-P'n 

+ . 2 i 247 
(n = 0,1,2,...) are [ef. (7), Appendix III] }a-tr-tnye~™*, (mr)-te-1'", 
(47)*"i"erfe(n/274) respectively, where 

ierfe x = erfex 

<n 


"| 


erfe x dx. 


and i"erfe x | 


Zz 


TABLE 1 


Values of Ay ms, from ref. (3) 








’ = 73 
20160 40320 
3 


161280 
































In this Table the constant A is given by A = ¢”(0) = 1-2326. 
Minor numerical errors in some of the coefficients from equations (25) to (33) of ref. 3 
are corrected: Mr. Glauert has kindly agreed with these corrections. 


If f(y, 7) denotes the inverse of /(7, p), then by formally inverting (4.8) term 
by term we obtain 


a ine) 
Sat) = te D Gam 7 ete verte Y Omnia )T He + 
m= ’ m=0 


i 2) w 
| i 

> 2" ya 

n=0 m=0 


)rininerfe(n/2r4) 


m 
mm+n+2 7) 


A,(n)r-te-7" 47 + A()r-te- 1°47 + s B,(n)rinierfe “+ (4.9) 
n=0 aT 


On inverting (4.5) we obtain ((7), Theorem VI) 


T 


F(y,7) = | F,(r—s)f(n, 8) ds. 


0 
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The skin friction is given by (3.4) as 


rata) (4.11) 


where 7,,, = “a'v-txg¢"(0), and from (4.10) with (@F/én),.9 written as a 
limit we have 


F fr) ————~ 


# — s)y/. l 41,2 G1 n+ a= 
a. (F,(7) F(r— 5 = + Gat > r+ j (4.12) 


0 


in which the coefficients a,, ,, are given in Table 1. 
Consider now the particular problem of a wall at rest, impulsively set 
in uniform motion at rt = 0 with velocity «, that is, F,(r) equals the unit 


function H(r). Then the velocity distribution becomes, by (4.10), (4.9), 


F(n,7) | f(n, 8) ds 


T T 


A,(n) [ ster ds-+ Ag(n) [ stems ds+ 
0 


0 


+ - Bn) | sininerfe( 5) ds 
n=0 0 


2atn-1A,(n) erfe A—z*yA,(n)(erfe A—7 te") + 


+ S (—)"2-"n"+2B, (n) —"—, (4.13) 


n=0 


wi 


where A = »/2rt and 
J, = (—)"4(n+2)! f w™-tinerfow du. 
When n > | it is found by repeated integration by parts that 
iu +4 S (—Ye+Dia+-tiverfo a. (4.14) 


Having determined Jj, J,..., J; from (4.14) and having evaluated J, we 
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find that 


= (1+ }A~*)erfe A—2-te-"A-1 

= (1+ 3A~)erfe A—2-te-"(A-1-++-A-3) 

Jy = (14+3A-2+ $A—)erfe A— 7 ¥e-"(A-14. 0-3) 
= (14+5A-2415\-4*)erfe A—a-he dy 149-84 2-5) (415) 
= (1+ $A-2+4-4A-4-+- 8-8 )erfic A— a —te-*(A-1-4. 70-3 4.885) 

Js, = (1+-34A-2-4.105)-44 195-6 )erfic A— 





—n-te-*(A-14. 10A-3-4- 87-54. 6-7) | 


Hence the first eight terms in the expansion (4.13) have been determined, 
A,, A,, and B, being defined by (4.9), the a,,,, being given in Table 1 
and the J,, being given by (4.15), A denoting »/2r#. 

The corresponding skin friction is given by (4.11) and (4.12) and the 
integral in (4,12) vanishes since F,,(r) and F,(r—s) are both unity through- 
out the range of integration. Hence 


x 


" 1+in 
(=) ee ae 2a {=}? + 7... (4.16) 
©1} n=0 (77)! “\ n=0 a+ 9) 


4.2. Solution for large times 


We now expand f(y, p) in a series for small p and invert term by term 
to obtain a formal solution for large 7 (Carslaw and Jaeger (7)). Glauert 
(3) has obtained a solution for small values of p in the form 


F(n.P) > P'S), (4.17) 
where fo(n) = $"(n)/¢"(0) 
7 © 
fur) = —4"| | PTS ardytT | bef, d»| 
0 


7 


= [ban (4.18) 
j 





” 
I = | e-#/6"2 dy 
F 


On inverting his solution (4.17) term by term, we obtain the result 





= d8(r) 
Sins) = > Sut) Pe (4.19) 
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where 8(r) is the delta function. Then F(n,7) is given by (4.10) so that 


F(n,7) TOS) as 











n 








= > ht | Ra(r—8) PO) ds = > fun) “32, (4.20) 


by successive integration by parts, and consequently the skin friction is 
given by equation (4.11) with 


(0 . 4.21 
(Fn) = > fi zee) (4.21) 


n 





In the particular case of a wall set impulsively in uniform motion from 
rest, for which F.(r) = H(r), the unit function, equations (4.20), (4.21) 


reduce to F(,7) = fol) | 


fal : 7 
pe = 0 
Fh, fi(0) 


for large +, which agrees with Glauert’s result for the wall in steady motion 
as we would expect. The value of f4(0) in (4.22) was given by Glauert, 


and is f,(0) = —0-8113. (4.23) 


(4.22) 


4.3. Summary of Laplace-transform solutions 

In sections 4.1 and 4.2, we have obtained a formal solution for the 
unsteady velocity distribution (3.4) which is given by (4.10), (4.9) for 
small times and by (4.18) and (4.20) for large times, «F,,(r) being the 
velocity of the wall. The skin friction is given by (4.13) and (4.15) for 
small times and by (4.11) and (4.21) for large times. 

In the particular case of a wall at rest, impulsively set in uniform motion 
at + = 0 with velocity e, the velocity for small times is determined from 
(4.13) to (4.15) while for large times it is determined from (4.22). The 
unsteady skin friction is proportional to —(@F/én),.9, which is plotted 
in Fig. 1. It is given by (4.16) for 7 < 1 and the limiting value as + + 00 
is given by (4.22) and (4.23). It is readily seen that for small times the 
curve very rapidly approaches the limiting value; so this curve may be 
extrapolated to its limiting value to a high degree of accuracy. Hence, 
for this particular flow, the skin friction has been obtained with sufficient 
accuracy by use of the Laplace-transform solutions only. In general, 
however, the curve for small times may not link up with the limiting value 
for large times. For this reason an alternative approximate method is 
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given in section 5, which may be used to achieve an optimum join between 
the solutions for small and large times. 

‘ satin 
Laplace transform method (§ 4) 


—-——-— Recommended curve for 7 >! 
- Pohihausen method (§5) 





— — 3 _ ert: 


r 


Fic. 1. Variation of unsteady part of skin friction with time. 


It is readily shown that (4.10) formally satisfies (3.5) provided that 
certain differential equations are satisfied by the functions A,(), A,(»), 
B,(m) of (4.9). It follows from these differential equations that, in place 
of the power series definitions, these functions may be written as functions 
of ¢ which are more amenable to calculation for values of » which are 
not small. We find 


T 


A,= 5 exP|—4 [ d ar] 
— 0 


78+ | i ¢? an)exp|—4 ( i$ dy] ’ 


St 
0 





7 7 
B, = &| —12¢’+762+3¢ Banta *an) | {4 (ddr) 
6 | fs n) exP|— | an} 


(4.24) 


and the rest are determined successively by the recurrence relation 


q 
By .exp{t [ $n} = —BA(0)+(B,+48B, expla { [4 an) — 
6 


7 


— J 9’ +1498, exp | {4 dn} dy. (4.25) 


In fact they are all of the form of the product of exp|—4 { dn} with 
fr) 
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a function of é which tends to infinity as » > 0 like a power of n. A rigorous 
justification that (4.10) satisfies (3.5) presents difficulties arising from the 
lack of knowledge of the properties of B,(n) for large n. However, it 
appears from the case of an impulse, that the Laplace-transform method 
for small times yields the exact solution. Similar comments can be made 
for the large times result, (4.20). 

For the impulse case we define a time, 7,, which is characteristic of the 
time taken for the flow to settle down into its new steady state and may 
be defined as follows. Let 7, be the steady limiting value of the unsteady 
skin friction. Then the characteristic time, 7,, is defined by 


t= | (‘t= —1) dr, (4.26) 
Ts 
0 
and this integral may be evaluated numerically to yield 
7, = 0-60(8). (4.27) 


The time taken for the flow to settle down to its ultimate state will be 
about three times this. Since + = at and the undisturbed displacement 


thickness, 5,9, is given by 
- 


(a/v)'845 = | (1—¢’) dn = 0-6479, 
0 
it follows that vt,/824 = vr,/adz, = 1-4(5). (4.28) 
It is noteworthy that this dimensionless characteristic time is of the same 
order of magnitude as the characteristic time for the asymptotic suction 
profile (Watson (5)). In that case vt,/5?, was unity. For a boundary layer 
of displacement thickness 0-1 cm, in air, t, is here about 0-1 sec. 


5. Alternative approximate method of solution 

We now obtain an approximate solution of (3.5), (3.6) by a Pohlhausen 
type of method. In this we satisfy the differential equation (3.5) only at 
the wall, at infinity, and on an average by integrating it across the 
boundary layer. When integrated with respect to 7 from 0 to 00, equation 


(3.5) becomes ad ~ 


[ Fan = Ur ry-0-2 [oF ay (5.1) 
0 0 


where use has been made of the boundary conditions (3.2) and (3.6). We 
now approximate to F by means of a cubic in 7/A, 


F = A+ B(n/A)+C(n/AP+D(n/A® forO0 <_< A 
and by | (5.2) 


F=0 for7 >A 
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where A, B, C, D are chosen to satisfy the four boundary conditions 
F(0,7)= F(t)  F(A,7r) = 0 \ 
dF,,(7) 
d 


> 


' 5.3) 
F’(0,7) = F’(A, 7) 0] ( 
These conditions lead to 

dF 


de n/A)—2(/A)?+(9/A)} 


F(», 7) $F {2—3(7 A)+(» A)j?}—}A?2 
(5.4) 
for 0 < » < A. We note that this is the typical single-parameter velocity 
distribution which arises in Pohlhausen treatments, the parameter, 
Fae 
FF...) v 
form. We now substitute (5.2) and (5.4) into (5.1) and obtain the ordinary 
differential equation for A, 


a(6r, AQ <\ 1A s & F.—6A dF, , AB ‘a Y 
( 


dr } dr = dr dr? 


du , = : 
—1, being similar in 


52 
? and the Pohlhausen parameter, A ( ‘ j 
at 


aa 


A 
i ‘ , oF, . , ; 
96 Bet? 3(n/A)+(y A)}—}A2° = {(n/A)—2(/A)? + (» silo dn, 
at 

0 

(5.5) 
which can be solved by standard methods. Having determined A as a 
function of +, we find from (5.4) that, for + > 0, 


eF\ i, (5 dF, (5.6) 
(: <4 0 2A 4) dr’ , 


a known function of 7. The skin friction follows from (4.13). 
We now give details for the particular case considered earlier, namely 
when the wall at rest is impulsively set in uniform motion at 7 = 0, so 


that F = 0 for 7 < 0 and F,(7) is constant for t > 0. Then, for + > 0, 
(5.5) reduces to 
1dA ae > 
+ a. — { 3) oh’ 
4dr ee {1 —3(/A)+4(n/A) } d' dy. 
0 
We may further approximate to the steady velocity distribution, 
Uo/Uy = $'(n), 
by its Pohlhausen value. Hence, following Goldstein ((8), p. 156), 
¢ = 1—(1—n*)3{1+(1—JA)n*} for n* <1) 
“=1 for n* > 1) 
n* = 2/8 = /Ab, A = (8/v)dU,/dx = 8*a/v 
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and for stagnation point flow, A = 7-05. We assume, as is reasonable 
physically, that the unsteady boundary-layer thickness, (v/a)!A, will not 
exceed the steady boundary-layer thickness, 5, that is, (v/a)!A <8 or 
A < Ai = (7-05)! for all values of +. This has been verified a posteriori 
because A is found to vary monotonically between 0 and 2-079. It follows 
that we may use the quartic expression (5.8) for ¢’ throughout the un- 
steady boundary layer. Then on substituting into (5.7) and integrating 
we obtain an ordinary differential equation which integrates immediately 
to give 


AdA 
#f1/, A\AP AP 3 (, AVAS 1/, AVS)’ 
30 +6 At 48 i40\" 3}ait gol! 6/ A? 
(5.10) 


the limits of integration being chosen so that A = 0 when 7 = 0. The 
integrand of (5.10) may be split into partial fractions and integrated 
exactly, yielding an analytic function + of A, or A may be found as a 


function of + from (5.10) by numerical integration. The quantity (=) 


then follows from (5.6) with F,, = 1, namely 





n=0 


é 3 
onan aD eee 5.11 
3 ee mn 


This is plotted in Fig. 1. 

We note that the curve for the unsteady skin friction as determined 
by the Pohlhausen method shows good agreement with those obtained by 
the Laplace-transform method for small and large times. In general the 
former will enable an interpolation to be made in any range of values of + 
for which the latter curves are not valid, though clearly in the present 
case satisfactory interpolation is possible without the use of the Pohlhausen 
method. 

It can readily be shown that, for small times, the unsteady skin friction 
as given by the Pohlhausen treatment has the correct form (in the sense 
that it varies like +~+) but is 6 per cent too low compared with the leading 


term of the exact expansion (4.16), whereas at large times the limiting 


value is 11 per cent too low. 3 
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SUMMARY 


The stream function for the shear flow with hyperbolic velocity profile past an 
elliptic cylinder has been determined as an infinite series of Mathieu functions. It 
is found that the stagnation streamline of the flow is displaced towards a region of 
higher velocity, this displacement increasing (i) with increase of angle of incidence 
of the cylinder to the main stream, (ii) as the stream becomes progressively non- 
uniform, (iii) with increase of minor axis length when the major axis length remains 
invariant. In each case the displacement reaches a limiting value as the cylinder 
moves away from the axis of symmetry of the stream. These limiting values are 
reached at critical distances from the axis of symmetry, which decrease as the stream 
becomes progressively non-uniform, but these distances are approximately in- 
dependent of incidence. 

The pressure coefficients and the resultant force and moment coefficients associated 


with the cylinder have also been obtained, and investigated numerically for the flat 
plate type of cylinder. 


1. Introduction 


WHEN the vorticity is constant throughout the fluid, the characteristics 
of the flow in the presence of cylindrical bodies of simple cross-sections 
have been investigated by Tsien (1), James (2), and Mitchell and Murray 
(3). A close investigation of the stagnation points in the field of flow was 
carried out in (1) and (3), and in particular in (3) the displacement of the 
stagnation streamline at large distances from the cylinder was determined. 
Considerable interest has for some time been concentrated on the pitot- 
tube effect concerning the displacement of the stagnation streamline by 
a pitot-tube set in a shear flow. Although this is strictly a three-dimen- 
sional effect as shown by Hall (4) and Lighthill (5), the two-dimensional 
approach of Murray and Mitchell (6) as applied to a circular cylinder 
representation of the pitot-tube confirms the main experimental result 
due to Young and Maas (7), viz. that the displacement of the streamline 
is towards a region of higher velocity. However, it is found that the two- 
dimensional theory underestimates the magnitude of the displacement as 
compared with experimental results. 

The flow of fluid past a thin aerofoil in a shear flow with parabolic 
velocity profile was previously investigated by the author (8) with the 
main object of determining the hydrodynamical forces acting on the 

(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 2, 1959) 





192 E. E. JONES 


aerofoil. The solution was an approximate one, the results being applicable 
only when the aerofoil was in the vicinity of the axis of symmetry of the 
field of flow. The present study is an extension which gives an exact 
treatment for a more general type of shear flow past an elliptic cylinder 
and a flat plate, and thus enlarges on that covered by (6). Since the 
subject-matter of this paper was completed a further paper by Murray (9) 
has appeared which investigates the non-uniform shear flow of a fluid 
past cylinders with sections of a general shape, the elliptic cylinder being 
a particular example. However, the analysis is complicated and the 
determination of certain coefficients tedious to carry out, even in the case 
considered when the cylinder has its axes parallel to or perpendicular to 
the direction of flow. It was noted, however, by Murray that an analysis 
involving Mathieu functions was possible, although this was not developed. 

In view of the fact that the analysis of this paper involves Mathieu 
functions many references are made to the treatise of McLachlan dealing 
with these functions, (10). The numerical analysis has been based on the 
N.B.S. tables (11), with reference also to Watson’s treatise on Bessel 
functions (12). 


2. The stream function 

The fluid is assumed to be incompressible and inviscid, and the flow is 
two-dimensional, steady and rotational, i.e. it is associated with a vorticity 
which in general varies with position in the fluid. Attention is confined 
to any one plane of flow and a set of orthogonal Cartesian axes O'X, O'Y 
locate position in this plane. The fluid velocity has components (vx, vy) 
relative to these stream axes, and the continuity equation can be satisfied 
by introducing a stream function 4(X, Y), where 

Vy = O/eY, Up = —Op/eX. 
The vorticity vector has its non-zero component directed normal to the 
plane and is of magnitude 
Up ov 


L = — - x — _ 2 ] 
~ a oe i ©) 


In steady motion the vorticity is constant along a streamline, implying 
V2 = f(b), where f() is an arbitrary function of %. Following the method 
of (6) and (9) a possible rotational flow is defined by choosing f(y) = w/c*, 
where c is an arbitrary standard of length, thus leading to the flow equation 


V2b = w/c. (2) 


In the absence of the cylinder the undisturbed stream is assumed to be 
defined by a stream function ¥,(X,Y), the vorticity in the stream being 
defined by (1) and the flow equation given by (2), with % replaced by yp. 
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It is assumed that the streamlines of the undisturbed flow are all parallel 
to the axis O’X, hence yf, is a function of Y only, and the flow equation 
can be written as d*s,/dY? = y,/c*, with a solution 


r 


he = —eV (sinh : +N cosh}, (3) 


Here V is the stream velocity along the negative X-axis, and N = cw,/V 
is a dimensionless constant, where w, is the vorticity component of the 
flow on the X-axis. If N = 0, then the flow is symmetrical about the 
X-axis, the velocity profile being of the hyperbolic cosine form, with 
velocity V along the negative X-axis and zero vorticity on this axis. The 
lengtheis a parameter indicating the deviation of the stream from the 
uniform. The stream is similar to that which occurs behind a symmetrical 
obstacle set in a uniform flow, and will be the main subject of the present 
investigation. 

Ultimately a cylinder is to be introduced into the field of flow with 
generators perpendicular to the flow plane. The right section of the 
cylinder has an elliptic profile, and the cylinder axes Ox, Oy are taken 
respectively along the major and minor axes of the ellipse with O at its 
centre. The stream function for the disturbed flow past the cylinder is 
Y(x,y), where v, = &b/ey, v, = —é@b/ex are the velocity components in 
the direction of the cylinder axes. As for the undisturbed stream it follows 
that (2) is satisfied by the stream function % in terms of z and y. If O 
has coordinates (Xy,¥,) relative to the stream axes, and the positive 2- 
axis is inclined at an angle @ to the positive X-axis, then 

Y = Y,+2sin@+ycos6. (4) 

When the cylinder is introduced in the stream the stream function yf, 
is modified by an amount V, becoming ¢ = 4%,+‘Y. Since both y, and 
ws satisfy a flow equation similar to (2), then 

VY = ¥/e*. (5) 
It is noticed that when ¥Y + 0, then V4 — V*f,. This verifies that at 
large distances from the cylinder where the latter has little effect the 
disturbed stream has the same hydrodynamical characteristics as the 
undisturbed stream. 

The elliptic profile of the cylinder section is assumed to have a semi- 
major axis of length u = hcosh € and a semi-minor axis of length h sinh 5, 
where 2h is the focal distance of the ellipse, and £, > 0. The natural 
coordinates for flow past the ellipse are the elliptic coordinates (£, »), where 


x = heosh€ cos », y = hsinh é sin , (6) 


with € = &, on the ellipse, and 0 < » < 2z. 
5092 .46 oO 





194 E. E. JONES 


With this change of coordinates, (5) transforms to 
rT .0F . 
-4+—— = 2q(cosh 2—cos 2n)¥, 
with q = k? = h?/4c?, and if the function is expressed in the form 
Y = fi(n)f.(€), then on separating the variables there results 


d*f, 


+ + (a+-2q cos 2n)f; 0, (7) 
dy? 
as -(a+-2q cosh 2£)f, = 9, (8) 
dé 


where a is the separation constant. These equations are respectively the 
ordinary and modified Mathieu differential equations (10). It is necessary 
that W(é,7) be periodic in » with period 27, and that it must tend to 
zero as € > +0, i.e. at large distances from the cylinder, where ys > ys. 
These conditions are satisfied by a suitable choice of f,(7) and f,(é). The 
differential equations (7) and (8) have solutions which satisfy the required 
conditions for a discrete set of eigenvalues of the constant a. 

The functions f(y) are then the ordinary Mathieu functions of integral 
order denoted by ce,,,(y, —q) and se,,(7, —q) with m > 0 (10, p. 21). These 
are continuous functions for all real », and their series representations are 
absolutely and uniformly convergent. 

Solutions for f,(€) take the form Fek,,(€,—q) and Gek,,(£, —q), which 
are the modified Mathieu functions of integral order (10, p. 248). These 
functions converge rapidly and monotonically to zero as > +00, the 
rate of convergence increasing with increase of £. The series representations 
converge uniformly in any finite part of the €-plane, including € = 0. 

The two sets of functions form solution pairs (10, p. 174), the complete 
solution of (5) in terms of elliptic coordinates being 


Ms 


WicV = { D,, Fek,,(€, —q)ce,,(n, —q)+ E,, Gek,,(€, —q)se,,(7, —9@)}, 


(9) 


m=0 


where D,, and E,, are dimensionless arbitrary constants, with EZ, = 0. 


3. Determination of the coefficients 

The coefficients D,, and E,, for m = 0, 1,..., can be determined by im- 
posing the boundary conditions at the surface of the cylinder, viz. that 
the profile of the cylinder section is a streamline of flow. Mathematically 
it is necessary that %/cV = p, a constant, when € = &, for all 7 in the 
range 0 < 7 < 27. Now k = h/2c, and if we write 


w = cosh £ cos n sin @+sinh é sin 7 cos 8, 
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then in terms of elliptic coordinates, with Y defined by (4), the undisturbed 
stream function y, takes the form 


jy = —eV{sinh(b+ 2kw)+N cosh(b+ 2kw)}, (10) 


where 6 = Y,/c. Ultimately it will be necessary to express quantities in 
terms of the ratio L = Y,/u, ie. L = 68, where 8 = c/u. All physical 
quantities for the flow past the cylinder can thus be expressed in terms 
of the non-dimensional quantities NV, q (= k*), 8, and L (or b). It is noted 
that € = sech~(2k8). The boundary condition on the cylinder requires 
that Yy+¥ = cVp, when & = &5, or from (9) and (10), 

p+sinh(b+ 2kw»)+ N cosh(b+ 2kw») 


> (Dn Fek,,,(&, —q)ce,,(7, —q)+ E,, Gek,, (&>, —q)se,,(7, —q)}, (11) 


m=0 

where w, = {w}¢_¢,, and is to be satisfied for all in the range 0 < n < 2z. 

It is possible to determine D,, and £,, as functions of p by using the 
normalized orthogonality relations for the ordinary Mathieu functions 
(10, p. 24), which are 
oe 7 (m=), 
| ce,,(7, —q)ce,(n, —q) dy — ‘ 

0 (m ~n); 


0 


[anda titade--oed q aiiteinn 
8€,,(7, — 7)8€,(7, —G) an = |. : 
} O (mn); 
2n 


( ce,,(7, —9)se,(n, —q) dyn = 0 (all m, n). 
0 


On applying these relations to (11), it can be shown that 
D,,, Feke,(€, —G) = 2(—1)"pA3"+-2(sinh b+ N cosh b)C,,,(k), 
Dyn FeK on 44(€o, —9) = 2(cosh b+-N sinh b)C,,, ,,(k), 
Ban sg GK gn, 42(€, —7) = 2(sinh b+ N cosh b)S,,, ,2(k), 
Bay. GK an 4(€, —7) = 2(cosh b+N sinh b)S,,, ,,(k), 


where 


2a 
C,(k) = (— 1)"C,,(—k) — ~ | en ce, (7, —q) dn, 
0 


2n 


S,(k) = (—1)"Sq(—b) = 5 | eth se,n(n, —9) dr. 


1. 
Expressions for these integrals are determined in the appendix. 
Equations (12)-(15) determine D,,, in terms of p, and Dj, 41, Lense, 
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E,,,.,, in terms of known functions. If p is known the function ‘’, given 
by (9), and thus ¢ = ¥,+‘Y, are known explicitly. The flow pattern is 
thus uniquely determined. 

The boundary condition already imposed on the stream function is not 
sufficient to determine the constant p. The ideal state of affairs would be 
to determine the lift force on the cylinder theoretically in terms of p and 
compare this with experimental values for the lift, thus determining a 
value or range of values for p. In order to obtain a theoretical result to 
the problem recourse can be made in the case of the elliptic cylinder to 
the assumption that the circulation round the curve coinciding with the 
section profile of the cylinder in the disturbed stream has the same value 
as that round the same curve in the undisturbed stream, cf. Tsien (1). 
It must be pointed out, however, that this extra condition, although a 
plausible assumption, has been introduced in order to provide a solution 
to the problem, and in fact requires experimental justification. 

If s is the curved profile of the cylinder section then the Tsien condition 
becomes , 

| (OV /e€)e_¢, dy = 0. (16) 


. 
8s 


From the definitions of the ordinary Mathieu functions, 
2a J—(_.])n 42n » * 
f ~~" {27(—1)"Ag” (am 2n), 
m| \0 (m = 2n+-1) 


0 
and | se,,(7,—q) dy = 0 (all m). 
0 


Hence substituting for ‘Y from (9) into (16), it follows that 


x 


> (—1)"D,,, Ag” Fek;,, (£5, —4) = 9, (17) 


m=0 
where (’) signifies differentiation with respect to £. The constant p is now 
determined by substituting for D,,, from (12) into (17), leading to the 
result 


p = —K(sinhb+N coshd), 
2, (— "Ao" Fan(bo — DC an(k) 
where K(é,, 9,q) = *=*— _—_—______—— (18) 
S (AB Fall —q) 


and by definition 


Fyn (€o; —q) - Fek; 2n (fo, — q) )/Fek,,,( Eo, —4q). 


Substitution for p into (12) gives the final form of D,,,. If this expression 
for D,,, and the coefficients in (13), (14), and (15) are substituted into (9) 








THE ELLIPTIC CYLINDER IN A SHEAR FLOW 197 


this yields the final expression for Y’, and with %, given by (3) the stream 
function ys of the flow past the elliptic cylinder. 

It is possible to approximate to the parameter K in (18), and thus to 
the stream function y, for streams with small velocity gradients. In this 
case the stream is nearly uniform, and then 8 = c/u is large, or q is small, 

The number q is limited such that 1 > \qlogq| > q. Then it is known 
that A! ~ O(g"-") (10, p. 46; Errata p. xii), and from (10, p. 382) it 
can be shown that F,,,(&,—q)/Fy(€,—¢q) is O(logq). Again from (10, 
p. 46) it is seen that A? = }q+ O(q*), Af = —4qA$+ O(q*), and using the 
series forms for the modified Mathieu and Bessel functions it can be proved 


-- FE, —q) __2(flogg-+y—log 2+) 


Fy(fo,—@) 1 + }q(log q)sinh 2&, 
where y is the Euler constant. 


On using the result J,,(z) ~ O(q"), it can readily be shown that (18) re- 
duces to 





+ O(q log q), 


ge) 
4S 

+-:4q*{(cosh 2£,—cos 20)? + 2(1 — 2+ 2£,)(1—cosh 2£, cos 20)} + 

+ O(¢ log q), (19) 


where » = log2—y = 0-1159..., and S = 1+-$q(logq)sinh 2£,. 
These results simplify in two limiting cases: 


K = 1+ -4¢(cosh 2€,—cos 26) + (cosech 2£,—coth 2€, cos 20) + 


(i) If h-> 0 and &, > © simultaneously such that he + 2R, then the 
ellipse degenerates into the circle of radius R (10, pp. 367-9). It can 
readily be shown from (18) that in general K — J,( R/c), a modified function 
of the first kind, and that % reduces to that deduced in (6). In particular 
for the near-uniform stream, K = 1+ R?/4c?+-O{( R/c)*}. 

(ii) If &, > 0 then the cylinder takes the form of a flat plate, and from 
(A 1) of the appendix, z = 2ksin@ and 6 = 32. It follows from (A 5) and 
(A 6) that S,,(k) = 0 for m > 1, and from (14) and (15) it is seen that EZ, = 0 
form > 1 in (9). For the flat plate of length 2u, with £ = 0, (19) becomes 


K = 1+ qsin®é-+ }q*(log q)sin?6 + }q?(1— 2u+-sin?@)sin?0 + O(q° log q). 


In this case of the flat plate the stream velocities at both leading and 
trailing edges become infinite in magnitude. This means that in the 
vicinity of these edges the stream function as deduced using the Tsien 
condition will not define a real fluid flow. 


4. Displacement of the stagnation streamline 
The displacement of the stagnation streamline from the cylinder for the 
general type of stream with hyperbolic velocity profile can be determined 
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by use of the stream function given by (9). It is then possible to determine 
the effect of cylinder shape, cylinder incidence, and form of the stream 
on the displacement of the stagnation streamline. This streamline curves 
off to infinity when it leaves the cylinder, its equation being &(£, 7) = cpV. 
At large distances from the cylinder ys —> yo, hence if the stagnation stream- 
line there has a displacement from the X-axis of magnitude d, its equation 
from (3) is given by 


— » d is . 
sinh —-+-.V cosh A (sinh 6+ N cosh bd). (20) 
c c 
On solving this equation we find 


(N+ 1)2 log| K (sinh b+N cosh b)+{K?(sinh b+ N cosh b)? + 1— .N2!4]. 
> 


Since the numerical analysis which follows refers to cylinders of invariant 
major-axis length, and the stream form is allowed to vary, it is more 
convenient to express the displacement of the stagnation streamline as 
a multiple of the semi-major axis length of the elliptic section of the 
cylinder, viz. u. If the displacement of the streamline is measured from 
a line through O, the centre of the cylinder, parallel to the X-axis, i.e. the 
line Y Y,, then the parameter determining its magnitude is 


D (d—Y,)/u, or D a(< —»), 


c 


where 8 = c/u, and dc is given by (21). 
When the stream is symmetrical about the X-axis, with V = 0, then 
from (21), ] 
D = Bsinh (K sinh 5] L, (22) 
B 
where L = Y,/u = fb. It is noticed from (22) that D(L) D(—L), 


hence at mirror image points in the X-axis the displacements are equal in 
magnitude but oppositely directed. Also from (18) it can be shown that 
K(@) K(2—6), hence from (22) it is seen that D(@) D(xz—8@), and it 
is thus necessary to consider incidence only in the range 0 < @ < }z. 
When the cylinder is at large distances from the X-axis, i.e. as L > co, 
then the streamline has a limiting displacement given by D, = Blog K. 
Also when the centre of the cylinder is on the X-axis, ie. L = 0, then 
T, = (dD/dL);-, = K—1. Other results of interest are (i) for given 8 
and L, dD/dK > 0 according as L > 0, i.e. at a given position of the 
cylinder away from Y = 0 in a given stream, D increases with K, (ii) if 
dK /d@ > 0, dD/d@ > 0 according as L > 0, indicating how the displace- 
ment varies with incidence at a given position in the stream, (iii) for 
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given 8 and K,dD/dL > 0 according as K > 1, alsodD/dL > 0as L>o, 
which indicates how the displacement varies with position of the cylinder 
in a given stream. 

In order to indicate the order of magnitude of D, the particular example 
of an elliptic cylinder in a symmetrical stream is considered. If the minor 
axis of the ellipse is half the length of the major axis, then tanh £ = 4 
and h = 4\3u. The stream parameter 8 = }v3, hence q = 1, and the 
stream has large shear. Values of K are determined from (18) with 6 = 0, 
kn, 4x. The functions F,,(&,—q) are derived from tables (11), and the 
C,,(k) as given in the appendix, the convergence being sufficiently rapid 
for computational purposes. It can be shown that 
K — 1,(z)+-0-15754 I,(z)cos 25 —0-00230 I,(z)cos 45 — 0-00006 J,(z)cos 68+... 


where tand = 2tan@, z* = 3(5—3 cos 26), and Jj, J,,... are modified Bessel 
functions of the first kind. Table 1 gives values indicating the variation 


TABLE 1 
Elliptic cylinder 





0 K D, T. 





° 1°39! 0°143 0-391 
‘7 1°972 0°294 o'972 
dn 2°689 0428 1689 














of K, D,, and T, with 6. The variation of D with L for the three stated 
angles of incidence is depicted in Fig. 1. 

For a circular cylinder with radius equal in length to the semi-major 
axis of the above elliptic cylinder, i.e. R = u, set in the same stream 
defined by 8 = }v3, it is found that K = 2-849, D, = 0-453, T, = 1-849, 
and these values can be compared with those recorded in Table 1. 

When the cylinder has the form of a flat plate, then €,—> 0, and 

2ksin@,& — 47. Table 2 gives the values of K, D,, and T, for angles 


TABLE 2 


Flat plate 


6 K Dt T. 








° 1°000 0°000 0°000 
hr 1°657 0219 0°657 
tn 2°530 0°402 1°530 














of incidence 6 = 0, }1, 47 assuming that 8 = }v3, and the length of the 
plate is 2u. Then since h = u, it follows that q = {in this case. In order 
to indicate the effect of cylinder shape on the displacement of the stagna- 
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Fic. 2. Cylinders at incidence @ = }2 with B 


tion streamline, in Fig. 2 the deflexion parameter D is plotted against the 
cylinder position parameter L for the flat plate, elliptic cylinder, and 


circular cylinder at incidence }7. 


2C 


) 


bv3. 
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For streams which are only slightly non-uniform it was shown in 
section 3 for small g, that K = 1+ O(q), hence for such streams from (22) 
it follows that the corresponding approximation to D becomes 


D = B(K— |) —(K—1 anh? tanh + Of(K— 1)%}. (23) 


This result has been investigated for streams with 8 = $v3 and 8 = 25v3, 
and it is found that change in incidence and cylinder shape follow the 





1B =V3/4N, 











8 = 5/3/25 









































6 ut 8 10 


Fic. 3. —— Elliptic cylinder at incidence 9 = }x. 


——-— Circular cylinder. 


same pattern as Figs. 1 and 2. The effect of change in stream form for 
a given cylinder at a given incidence is exhibited in Fig. 3—the cylinder 
has the elliptic section defined in conjunction with Fig. 1 and is at an 
incidence @ = }z. For comparison the effect of change in stream form on 
the displacement of the stagnation streamline from a circular cylinder of 
radius equal to the semi-major axis of the section of the elliptic cylinder 
has been included—this result is as given by the analysis of Murray and 
Mitchell (6), with modifications to allow a variation in stream form. 


5. The pressure coefficient 


The pressure equation for steady rotational flow in the absence of body 
forces (13, p. 109) is 


“+4t— [ V¥dp = 0, 
p 


AB 








202 E. E. JONES 


where C is an absolute constant, B is the point in the fluid at which the 
pressure P, the density p, and the velocity v are measured, and A is any 
fixed reference point. For the flow with hyperbolic velocity profile it is 
known that V2 = %/c?, and so 


5 (ye? pb?) C, (24) 


where ys is the value of the stream function at B. At large distances from 
the cylinder, i.e. where |X| -> co, the flow has the same characteristics as 
the undisturbed stream. If the pressure at such points is I], then by use 
of (3), it follows that 


p re 2c* 


The pressure coefficient associated with the cylinder is defined by the ratio 
2(P—I1)/pv%, where P is the pressure at points on the cylinder surface, 
and vy is the undisturbed stream velocity at the centre of the cylinder, 
derived from (3), with a value v, V (cosh b+.N sinh b), where b = Y,/c. 
At points on the cylinder 4 = cpV, where p K (sinh b+N cosh b), 


; N-+tanhb \? »\2 
C, = 1+(K? 1 : =) ©) ; (25) 
1+ tanh b Uo) 


By comparison the pressure coefficient C,,) for the cylinder in a uniform 


whence 


stream of undisturbed velocity v, is given by 1—(v/v,)*. 

We now determine v, the velocity in the disturbed stream on the surface 
of the cylinder. This is the tangential fluid velocity component with a 
value H~ab/e€, for the value € = &, and where 


2H? = h*(cosh 2£,—cos 27). 


By use of (10) it is possible to develop éy,/e&, when € = &, as a Fourier 
series with 7 as the variable, and in this connexion it is found expedient 
to use equations such as (A2) of the appendix. Similarly by using the 
series forms for ce,,(y,—q) and se,,(y,—q) in (9) it is also possible to 
develop é'¥/éé, when € = &, as a Fourier series with » as the variable. 


In fact 
l 


Cys) ~ . : 

: ( z) = Ay > (a, cosnyn +b, sinnn), (26) 
Vo CS é £o n=1 

where the a, and b,, are known constants depending on the geometry of 

the stream and the cylinder. It is comparatively easy by use of the above 
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methods to show for the flat plate in the symmetrical stream that 


a, = “tanh b> (1m 1)" KAR" —Cop AB Fam(0, —G) (n> O), 
m=O 
~~ : 
k P3 ( i pypremene 2m+1 


m=0 


Beet Femi, —G) (nm > 9), 


) 
= =" I,,(2ksin B)eot @tanhb (n > 1), 


2n+1 
bs, + P 


On introducing the constants Ay == dy, A 
is possible to deduce from (26) that 


I,,,.,(2k sin @)cot@ (nm > 0). 


- $(a,—ib,) for n > 1, it 


 ¥y 


I ey , ‘ inn7 —in 9 
73-37 = Hot J (une ™+fine'™), (27) 
h*v5 0€]¢ fo n=1 

where, for n > 0, 


n as 

oe y » 

Ln = > 4-An rt* 2, Be bgse 
r=0 r=1 


The result (27) can thus be used in (25) in order to determine C,,, since 


(29 a ceases an 
vy) (cosh 2£,—cos 2) \hvy O€) be . 





For the elliptic cylinder in a symmetrical stream it can be shown that 
eof] —cos 2(0+7)} | ktanh b 
cosh 2£,—cos2n cosh 2£,—cos 2 


eof 2 sin(@+ »)—sin 3(0+ )}—e&{sin(@+ 37) + sin(@—»)}]+ O(k*). 





[e-* sin(@+-»)— 
U] 


In order to avoid lengthy expressions this result has been developed 
only as far as terms of order O(k), and the result can be used to determine 
the effect of small stream shear on the pressure coefficient. The terms 
independent of k on the right-hand side of the equation give the value 
of C,,9, the pressure coefficient for the cylinder in a uniform stream (13, 
p. 161). The corresponding result for the flat plate is obtained on substi- 
tuting £ = 0. 

In order to indicate the order of magnitude of the change in pressure 
coefficient due to non-uniformity in the stream as well as change in 
incidence, the results for the flat plate in a symmetrical stream of large 
shear have been investigated—it is assumed that 8 = }v3, and hence that 
q = k* = 1. The coefficients tabulated at various values of the eccentric 
angle » for the angles of incidence 6 = }m and 4m are C,, C,,, and Cy», 
the pressure coefficients associated with the uniform stream, and the non- 
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uniform stream when L = 0-25 and | respectively. The values of C,, for 
larger values of L do not differ appreciably from those of C,,.. The infinite 
pressures at 7 0, 7, and 27 occur as expected, due to the presence of 
infinite fluid velocities at the trailing and leading edges of the plate. 


TABLE 3 
Pressure coefficients for the flat plate 














t 
0=4r 0=42 
U] C no Cn Cpe Cpo Cm Cpe 
x 23 x x x ° 8) 
7 I o°g2 1°33 ° 2°97 1°36 
4a o'5 1°43 2°30 I 1°33 18 
nr I 1°04 1°40 ° 4°45 1°43 
7 L x a > e) x oC 
q7 I I°21 2°6 ° 1°45 1°43 
$n o's 095 1°39 I 1°33 2°18 
jr I 1°31 I'go ° 2°07 1°36 
Tr x x x 20 x X 














The values for slightly non-uniform streams, when q is small, can 
similarly be obtained, but the pressure differences in these cases are not 
so pronounced. 


6. The force and moment coefficients 

(i) If the lift force on the cylinder in the direction normal to the main 
stream away from the X-axis is F;-, and the force in the direction of the 
negative X-axis is Fy, then 

F,,—iF, = e- [ P di. 
S 

Here @ is the angle of incidence of the cylinder to the main stream, P is 
the pressure on the cylinder of section contour S, and 7 = x—iy, where 
(x,y) are position coordinates referred to the cylinder axes Ox, Oy. The 
pressure P is defined by (24), and since % is constant on S, then 


Fy —iFy ke? | v* dz. 
s 


It is now possible to introduce force coefticients defined by 
Y . 9 ’ , 9 
Cy = Fy /zpur3, Cy = Fy/mpurj. 
From (6) it is seen that z = hceosh(€,+in) on S, and (v/v,)* is given by 


(28), hence 


27 


Cy—iC'y ] y cembihit. +-in)de. 29 
peers —— az | ( : 7) cosech(£,+-i) dy (29) 
0 
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But for €, > 0 it is possible to write 


we 
cosech(€,+i7) = 2 } e@n+iXfo+in), 
n=0 


and with the use of (27), integration of (29) leads to the result 
iCy = 2ie-" sech €, 2, Hon+1 e-em the, (30) 
n= 


For symmetrical streams with small shear, i.e. with q small, it can be 

shown for the elliptic cylinder that 

CO; = k(e%*—cos 20)sech &, tanh b+ O(k*), 
and the force coefticient Cy is at least of order O(k?). For the uniform 
stream it is seen that the force on the cylinder is zero. 

For the flat plate in a symmetrical stream higher order approximations 
can readily be determined under the limiting process €,— 0, and are 
recorded here as the expressions involved are not as long as those for the 
elliptic cylinder. In fact 

c,. ~ 2k{t 3 (4 sin29—1)\sin2@ Slog k 
; BT cs heen FAY —1)})sin?6é tanh b+ O(log k), 
A," 8 J 
where A, = 14+ (log k—p—4)+ O(k' log k), 
and Cy is at least of order O(k* log k). 

A numerical analysis of the results of (30) for a flat plate at angles of 
incidence 6 = 0, }7, 4a respectively in symmetrical streams of varying 
shear defined by 8 = 25v3, §v3, and }v3 respectively, has been carried out 
and values of Cy tanh 6 recorded in Table 4. 


TABLE 4 
Force coefficient Cy/tanh b for flat plate 


0 B = 25/3 $3 V3 








° ° ° ° 
jn o-o12 0120 1°634 
ad 0°023 0°242 18-256 














For constant f, i.e. for a given stream, the lift coefficient is proportional 
to tanh(L/), the constant of proportionality increasing with increase of 
angle of incidence in the interval (0,47). The relationship Cy—L for 
various values of 8 when the plate is at incidence }7 is exhibited in Fig. 4. 

At the other extreme of geometrical shape for the cylinder section, viz. 
the circular section, by use of the stream function of (6), the lift force 
can be expressed in the form 


Cy = 4Btanhb ¥ 1,(8-)/K,(B, 
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results for Cy-/tanhb are 0-046, 0-496, and 51-904. The force coefficient 
Cy is zero in magnitude. 


where 8 = c/R. For the three stream forms of Table 4 the corresponding 
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Fic. 4. Flat plate at incidence @ = }r. 


(ii) The couple acting on the cylinder in an anticlockwise direction has 
a moment M about the axis of the cylinder through O given by 


M = re Pz di. 


By use of (24), and the result 
re ( zdz=} ( d(zz) = 0, 
cS S 


it follows that the moment of the couple reduces to 
M re —tp ( vz a}. 
s 


If the moment coefficient is defined as C,, = M/mpu*v2, then as in (i) of 


this section 
2r 


; a ire 
Cy = re 5, sech*é, | ( 


2 
in) ¢ 3 
= “) coth(€,+in) dy. (31) 
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But for &, > 0, it is known that 


coth(£,+in) = 14+2 ¥ e-néorin), 
n=1 
hence on using (27), and integrating (31), there results 


a 
Cy = re 2isech’€, >, Han ents (32 
n= 


it being known that py is a real constant. 
For the elliptic cylinder in a symmetrical stream of small shear it can 
be shown that (32) reduces to the result (13, p. 165), 


Cy = }sin 20 sech*é,+ O(k*), 


and for the flat plate, allowing &, > 0, 


( 2 4 
-} 1+ <(16 sin?6— 3) + = [37-382 sin?6 + 576 sin*#} sin 26+- 


+ 4{k?+-}44(3 sin? + 1)}sin?@ sin 26 tanh*b+ O(k*), 
where 


Ay = 14 B(log k—p—})+ Blog k—p— J) + O(K* log k). 


A similar type of numerical analysis as applied in (i) of this section 
can also be applied here. It is found that for angles of incidence 6 = 0, 
kor the moment coefficient is zero in magnitude, and for @ = 47 the values 
of C\, are given by Table 5. 


TABLE 5 
Moment coefficient Cy, for flat plate 
B Cu 
25y3 0°500+0°000 tanh*) 


§V3 0521 —0°003 tanh*b 
4/3 0440+ 0°786 tanh*d 














The action of the couple is to orientate the plate to a position transverse 
to the main stream. The couple associated with a circular cylinder has 
zero moment about its axis. 


7. Conclusions 

Figs. 1, 2, and 3 indicate that as the cylinder is stationed further from 
the axis of symmetry of the stream, the non-dimensional displacement of 
the stagnation streamline increases from zero to a maximum limiting 
value, the latter being reached at critical distances which are approxi- 
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mately independent of incidence and shape of cylinder section but increase 
the more uniform the stream. For a fixed position of the cylinder the 
displacement increases (i) with incidence, the rate of increase being greater 
at smaller angles of incidence in the range 0 < @ < 4}n, (ii) the greater 
the transverse dimensions of the cylinder section, (iii) the more non- 
uniform the stream. In all cases the displacement of the stagnation 
streamline is towards a region of higher stream velocity. 

Table 3 illustrates how a change in stream form and change in incidence 
affects the pressure coefficient at various points on the flat plate type of 
cylinder. The effects are much more pronounced when the non-uniformity 
of the stream is more pronounced, i.e. for streams of greater shear, and 
the pressure distribution becomes asymmetric with respect to the centre 
point of the straight line section of the plate. 

Although by virtue of the complexity of the equations involved in the 
analysis it has not been found possible to prove that the force on the 
cylinder in the direction of the main stream is zero, it is conjectured that 
this is so in view of the results for the elliptic cylinder in a stream with 
small shear, and the numerical results for a flat plate in various streams 
at various angles of incidence. 

When the cylinder is at incidence @ = }7 to the main stream, it is acted 
on by a force in the direction of the major axis of its right section. A similar 
state of affairs exists when the stream has constant vorticity (2), and as 
noted by James (2) it is evident that as regards this transverse force effect 
the shear flow with distributed vorticity is analogous to the potential 
flow with circulation round the cylinder. 

A comparison of Figs. 3 and 4, and Table 5, shows that the lift force 
and moment coefficients exhibit the same general characteristics as the 
displacement of the stagnation streamline. The lift force coefticient for 
a symmetrical stream is directly proportional to tanhb, where b = Y,/c, 
and has a maximum limiting value which is greater for larger angles of 
incidence in the range 0 < @ < $7. It is noticed from Fig. 4 that the effect 
of an increase in the non-uniformity of the stream is to increase the 
magnitude of the lift coefficient. By comparison of the lift force on the 
flat plate and the circular cylinder, it is seen that the effect of increase 
of length of the semi-minor axis of the cylinder section is to increase 
the magnitude of the lift force at given incidence and position of the 
cylinder. 

The moment coefficient is a linear function of tanh’), and this also has 
a maximum limiting value. This coefficient is zero however at the ex- 
tremes of the range 0 < @ < 4m. There is an increase in the magnitude 
of the moment coefficient as the stream becomes more non-uniform. 
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APPENDIX 
In order to evaluate the integral C,,,(k) of (12) subsidiary variables z, 5 are intro- 
duced such that 
2k cosh €,sin@ = zsind, 2k sinh €,cos@ = zcosd. 


Hence tan 6 = cothé, tan@ 


2% = 29(cosh 2£,—cos 26) \’ 
and thus 2kw, = zsin(y +68) in C,,(k). But it is known that 


t sin(n +8) I,(z)+2 ¥ (—1)'I,,(z)eos 2r(n+ 5)+ 
1 


rT 


+2 2, (— War a(2)sin(2r +1)(n+8), (A2) 


hence on substituting the series representation for ces,(7, —9), 
" x 27 
C2,(k) on > (—1)"*"A? | e*sinin+Scos 2rn dn 


r=0 0 
(—1)" s A2",,(z)cos 2r8, (A 3) 
r=0 


where z and 6 are as defined in (A 1). 

If k changes sign then from (A 1) the angle 6 is replaced by 7+ in the above 
result, hence C,,,(k) = C,,(—k). 

In a similar manner it can be shown that 


© 
Conyilk) = (—1)" 2 Bot er41(z)sin(2r+ 1)8, (A4) 
r=0 


Sonuilk) = (—1)" & AP TTa,41(2)e0s(2r + 18, (A5) 


0 
Sonsalk) = (—1)" 5 BB F2I,,.9(z)sin(2r-+2)8, (A6) 
r=0 


the results for change in sign of k following immediately. 
Following an integral equation approach (10, p. 202) alternative expressions for 
(A 3)-(A 6) are determined in the form 


PanCan(k) = Cegn(Eo, —G)ceen(8, 9), 
PansiSangilk) = Seon41(Eo. —9)¢@an41(9, 9), 
S2n41Congilk) = Ceon,1(Eor —9)8¢ansi(9, 9), 
Seny2Sanyalk) = —Seon,2(For —9)8@2n42(9, 9), 
where the functions appearing in these results are defined in (10). 
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ON THE CALCULATION OF EDDY VISCOSITY AND 
HEAT TRANSFER IN A TURBULENT BOUNDARY 
LAYER NEAR A RAPIDLY ROTATING DISK 
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SUMMARY 

In this paper the distribution of the radial component of Reynolds shearing stress 
is calculated in the turbulent boundary layer over a rapidly rotating disk, by an 
integration of the equations of mean flow, the validity of the von Karman velocity 
profiles being assumed. The distribution of eddy heat diffusivity is then evaluated 
by applying Reynolds analogy in a thin region of flow very near the surface of the 
disk, and using the approximate similarity of the radial component of disk flow and 
flat plate flow to extend the results through the remainder of the inner, crucial, part 
of the boundary layer (about !6 per cent of the layer thickness). The method of 
calculation developed by Davies and Bourne (5) is then applied to evaluate, in the 
case of air, the rate of heat transfer from the disk, when the surface of the plate is 
kept at a constant temperature. It is assumed that the disk is rotating about its 
centre sufficiently rapidly to neglect (a) the completely laminar zone of flow near 
the centre, and (6) the thickness of the laminar sublayer; under these conditions the 
calculated value of heat transfer is found to be in good agreement with the value 
suggested by the experiments of Cobb and Saunders (4). 


1. Introduction 

Tue flow due to a disk rotating in its own plane was considered first by 
von Karman (1). He considered both the laminar and the turbulent cases 
and obtained, by his integral method, expressions for the velocity com- 
ponents in the boundary layer flow over the disk. The distributions of 
circumferential and radial velocity components have been measured experi- 
mentally by Gregory, Stuart, and Walker (2), who found agreement with 
von Karman’s distributions when the flow was laminar. Good agreement 
was obtained in the turbulent conditions for the circumferential velocity, 
and reasonably good agreement was also obtained for the radial component 
in the inner part of the boundary layer. 

In the laminar case the associated problem of calculating the rate of 
heat transfer has been solved by Millsaps and Pohlhausen (3), when the 
surface temperature of the disk is independent of position. No solution 
has been given, however, in the important case of turbulent flow, which 
occurs in most practical applications. 

Experimental results involving heat transfer into air have been obtained 
by Cobb and Saunders (4), when the flow over the disk was in part laminar 
and in part turbulent. It was shown that, as the Reynolds number in- 

(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 2, 1959) 
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creases, the transition point from laminar to turbulent flow moves inwards 
and an increasing proportion of the disk surface comes under turbulent 
conditions. Cobb and Saunders suggested that in the limiting case of a 
very rapidly spinning disk, when the completely laminar zone near the 
centre of the disk can be neglected, the experimental results approach 
values given by the expression 
N = 0-015 R°8, 

where N and R are appropriate Nusselt and Reynolds numbers for disk 
flow; here R = wr?/v, r and w denoting the radius and angular velocity 
of the disk, v the kinematic viscosity of the ambient fluid; 

N = Hr/k(T,—T), 
k denoting the thermal conductivity of the fluid, 7, and 7, the temperature 
of the disk and the ambient fluid, and H the average rate of heat transfer 
per unit area of the disk. 

In order to formulate a corresponding theoretical approach the analysis 
developed by Davies and Bourne (5), for heat transfer from a flat plate, 
is applied in this paper to the problem of a rapidly rotating disk. It is 
assumed that the disk is rotating sufficiently rapidly to neglect (a) the 
completely laminar zone of flow near the centre, and (b) the thickness of 
the laminar sublayer, these being the limiting conditions in which Kar- 
man’s velocity profiles apply. Reynolds stresses are then considerably 
greater than viscous stresses over the whole flow field. 

The equation involving the Reynolds stress component, associated with 
radial and normal eddy velocities, is first integrated, following the usual 
boundary layer approximations and applying the expressions given by 
von Karman for the radial and circumferential velocity components. The 
experimental results of Gregory, Stuart, and Walker (2) show that the 
position of the maximum radial component of velocity is given accurately 
by Karman’s expression and fairly good agreement is obtained in the 
inner part of the boundary layer. A considerable degree of error is found 
in the outer part of the layer. It should be noted, however, that a negligible 
completely laminar zone is assumed by von Karman, whereas the Rey- 
nolds number involved in the measurements is not sufficiently high to 
permit the completely laminar zone at the centre to be neglected. We 
note that near the surface von Karman assumes a one-seventh power law 
profile for the radial component of flow, because this is approximately 
valid in a turbulent boundary layer flow over a flat plate, and it is reason- 
able to suppose that the two flows are similar in the inner part of the 
boundary layers. The ensuing calculated values of Reynolds stress are 
expected to be accurate only in the inner part of the boundary layer 
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(about 16 per cent of the layer thickness), but, as shown by Davies and 
Bourne (5), this is the crucial region for heat transfer calculations. 

Using the concept of Reynolds analogy, which has been shown to be 
valid in forced convection from a flat plate, the distribution of eddy heat 
diffusivity, €,,, is then evaluated in the region of flow which is immediately 
adjacent to the surface of the disk (about 2 per cent of the thickness of 
the boundary layer). A power law representation of ¢,, in this inner flow 
is then constructed with good accuracy, and the dependence of this power 
law on distance from the surface of the disk is found to be in agreement 
with that given by a corresponding flat plate flow. Reynolds analogy is 
shown to break down at points further distant from the disk. However, 
since the von Karman profile of radial velocity for the disk flow is in close 
agreement with the flat plate profile in the inner 16 per cent or so of the 
boundary layer thickness, the power law for €, is assumed to be valid 
throughout this part (the inner 16 per cent) of the layer. 

The method of sources is next applied, the heated disk being regarded 
as an assembly of concentric circular sources. The radial distribution of 
source strength (i.e. of local rate of heat transfer) is determined by solving 
an integral equation, the temperature of the disk being taken to be constant 
over its surface. 

The calculated results are found to satisfy the relation 

N = 0-014R%, 


This is in good agreement with the result suggested by Cobb and Saunders 
(4). 


2. The equations of fully turbulent mean flow and of mean tempera- 

ture over a rotating disk 

In order to set up the Reynolds equations of mean flow in the fully 
turbulent case we neglect the viscosity terms in the laminar flow equations 
and, applying the usual Reynolds averaging technique, we obtain the 
equations of mean motion (see, for example, Goldstein, 6). When ex- 
pressed in cylindrical coordinates r (the distance from the centre of the 
disk) and z (the normal distance from the disk surface), these are 

y2 l (v'v’) 


gee we eee, a 
cr C2 r ror cz r 


a ev UV . Oe Woe” a 
pe ke iia in aw. pe 
U cr , cz + r > pre” ) aa” wv) r’ 


pel ye 


cr Cz cr 


a—, 8§—— uw 
= —2ww)-2 wv) -"*); 
} Cz r 
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U,V, and W denote mean velocities in the radial, circumferential, and 
normal directions respectively; u’, v’, w’ denote the corresponding com 
ponents of eddy velocity ; p(u'u’), p(v'ev’), p(w’), p(u'w), piu’), p(v'w ) 
denote the components of Reynolds stress, where p is the fluid density. 
The equation of continuity of the mean flow is 
l < : ow 
(rl) 0, (4) 
rol C2 
and the equation of mean temperature (neglecting friction heating and 


molecular heat transfer) is 
U W (u’' 7") (w’ 7”) 4 , (5) 


where 7” denotes the departure from the mean temperature, pc,(u' 7"), 
pc, (w 7") denote mean flux of heat in the radial and normal directions, 
and ¢,, is the specific heat of the fluid at constant pressure. In a turbulent 
boundary layer on a flat plate it is known from experiment (see, for 
example, Townsend, 7) that the various components of Reynolds stress 
are of the same order of magnitude, and, since the flow in the inner part 
of the boundary layer over the disk is similar to the flow over a flat plate, 
this result is probably valid in the inner part of disk flow. If we then 
make the usual boundary layer approximation and retain only terms on 
the right-hand sides of equations (1), (2), (3), and (5) involving derivatives 
with respect to z, we obtain equations involving only the three components 
of Reynolds stress associated with w’. 

The reduced forms of (1) and (5) only are required in the present calcula 


tion, and we obtain 


ol ,oU VY? . 
l + W “ (wu’), (6) 
cr C2 r C2 
and US d + Ww I : (u’ T”). (7) 
cr C2 Cz 


In equation (6), /’ and V are given by the von Karman profiles (1), viz. 


U rurE"7(1 —€), (8) 
and V = wr(l—€""), (9) 
where a = 0-526 and € = 2/4, 5 being the thickness of the boundary laver, 


as defined in Karman’s integral method of calculation, and given by 


ra) (v/a) V5p3/5, (10) 
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An expression for the vertical mean velocity component is now obtained 
from (4), which can be written in the form 
1eW eW Uo el 
5 of dz or Or’ 
where the right-hand side is a function of € only. By integration we 
obtain the form 


W wd 87 esr +5") (11) 
40° ri 
3. Calculation of the distribution of radial component of Reynolds 
shearing stress over a rotating disk 
Substitution of (8), (9), and (11) into equation (6) for the radial com- 
ponent of Reynolds stress, p(u’w’), and integration yields 


ae > l 7 . 
(u'w’) Pr _ y3¢4)9/5,,1/5,-8/5 _ £4. “£817 +. 
p . i 
ys ad ot agua 203 93/7 7\. (12) 
9\ 280 af T 600 725° J 


where +, is the radial component of shearing stress at the surface, and is 


"f 


derived by von Karman’s method (1) to be 

T, = 0-0225p(arw)7/4yV45-U4( 1 + a?) 98, (13) 
Tn a . . 
The caleulated values of (u'w’) in the inner boundary layer are shown in 
Fig. 1, demonstrating a change in sign at € = 0-042, from the zone domi- 
nated by the constraining effect of the surface of the disk in the radial 
direction to the zone dominated by the centrifugal effect. 


4. Calculation of heat transfer from a rapidly rotating disk 

In order to obtain a solution of the associated heat transfer problem 
the Reynolds analogy is now assumed to be valid in the region of flow 
immediately adjacent to the surface of the disk. We express this mathe- 
matically by the relations 


(uw ) (1 “w T”w’) 
——- == —=€, = —€p, 
eu /éz eT Jez - 
and € = €y,; (14) 


where ¢ and e,, are respectively the eddy viscosity and eddy heat diffusivity. 
Using (14), the mean temperature equation (7) becomes 
OT wT af8@ 
os wit = Stel in 
dz «= Oz\ Gz 


The left-hand side of equation (15) is then simplified by an appropriate 
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von Mises transform (8, p. 126), so that r and z are replaced by r and & 
(the stream function of the mean flow) as independent variables. We 
obtain 





a = #2 (uZ), (16) 
cr Cus Cus 

+0-004, 

+0002} 
a 
Re ° 

ls 

0-002 

-0:004 oe 4 — 
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. + + ae 
Fic. 1. Distribution of radial component of Reynolds stress, p(u'w’), 
near a very rapidly spinning disk. 


In order to apply the convenient method of calculating heat transfer em- 
ployed by Davies and Bourne (5), the term «/’ must now be expressed in 
terms of ys. An expression for % in terms of r and & is first obtained by 
integrating the relations 


—=rU and 


ous . Cds 
C2 cr 


TW, (17) 
which replace the continuity equation. Using the von Karman form (8) 
for LU’, this integration yields 

bb = a2w4/5yV5p13/5( 7 ¢8/7 _ Z £157), (18) 


An expression for «U in terms of r and é is then obtained from (8), (12), 
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(13), and (14): this is 
el = Ta®w5p25¢(1 —£)(1—8£)-rl/ T 0225a-"4(1 + a?)9/8+- 


| | 8/7 189 {7 20 16/7 203 23/7 iT 

io . i) 
‘| aft 4a? GE “9280 2 a)e"+ +oae T7358 (19) 
The variation of 8° with «U/r!¥5, as given by (18) and (19), is shown 
in Fig. 2 up to € = 0-04, 
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Fic. 2. Variation of eg U with J; y = eg U/(w8 vr), X = fb/(wi/yV5r9/5); 
caleulated from Reynolds analogy, ——— y = 0-0075X°, 





We note at this stage that all the available experimental evidence 
concerning forced convection in a turbulent boundary layer on a flat plate 
supports the Reynolds analogy method of linking through relation (14), 
momentum and heat transfer. It is found in the case of the disk flow that 
the calculated distribution of Reynolds stress in the thin region £ < 0-02 
leads to a distribution of « which is similar to that in flat plate flow. 
Consequently we have taken ¢« and e€,, to be identical for disk flow in this 
region. However, at € = 0-042, the calculated value of « is already zero, 
and the associated zero value then given by Reynolds analogy for €;, is 
not physically reasonable. There must be a mean, outward, non-zero flux 
of heat from the disk throughout the boundary layer together with a finite 
gradient of temperature at all points in the layer. Consequently relation 
(14) cannot be satisfied with «,, = 0, and Reynolds analogy breaks down 
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in this neighbourhood. Shear tlow turbulence theory, in its present stage 
of development, cannot be used to predict the form of €,, over the whole 
thickness of the boundary layer, and we make the assumption that the 
values of Ue,, calculated in the thin layer immediately adjacent to 
the disk surface can be extrapolated smoothly, through the part of the 
boundary layer which is likely to influence strongly the rate of heat transfer 
from the surface. This can be justified in the following way. 


We first suppose that e,, and uy can be related by a good power law 
PI H ; } 


representation in the inner part of the boundary layer £ < 0-02. We write 
y = O25, (20) 
where y € U/(w®5y25y1/5) and =X us /(qo8/ 5p/5p13/5)_ 


The variation of y with X is shown in Fig. 2: a very good representation, 
with b = 0-0075 and ¢ = 0-54 is found possible with 4 per cent deviation 
at most from the exact values in the region € < 0-02. As in the flat plate 
case (5), the ensuing calculated rates of heat transfer are found not to 
depend critically on the choice of 6 and ¢; i.e. any pair of values of 6 and f¢ 
which yields a good fit to the actual variation of y with X leads to caleu- 
lated rates of heat transfer in reasonable agreement with each other. The 
particular value ¢ 0-54 (with an appropriate 6 value) is taken, since it 
is also suggested by the following alternative approach. 

In considering the extension of the values of €n beyond é 0-02. we 
note that in the case of turbulent flow over a flat plate Davies and Bourne 
have shown (5) that a good power law representation of the €,, variation 
in the inner 16 per cent or so of the boundary layer thickness leads to 
calculated values of heat transfer in good agreement with experimental 
values. We infer that the flow properties in the outer region of the 
boundary layer do not influence significantly the rate of heat transfer at 
the surface; this is due to the fact that the mean velocity increases very 
sharply in the part of a turbulent boundary layer near the surface and 
consequently the main flux of heat takes place in this region. It is now 
reasonable to suppose that a similar result applies in the case of a rotating 
disk. In this problem we find that a one-seventh power law velocity profile 
constitutes a good approximaticn to the von Karman profile for radial 
flow in the inner 16 per cent or so of the boundary layer. It can be shown 
by calculation that the factor (1—€) in equation (8) can be effectively 
neglected in this part of the laver and (8) can be written in the form 

U U, £1/7 [z 1a(y ey )U 53/5} [1/7 (21) 
where wr, so that U/U, = 1 at € = 1. Comparing now (21) with 
the flat plate profile discussed by Davies and Bourne (5) 

U/U, = 9", 
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where 7 y (kx®) and x and y are downstream and normal coordinates, 
we see that the profile parameter denoted by qg in their work (5, equation 
(3)) has the corresponding value 3 in the disk problem. We then assume 
that the distributions of €,, in the radial component of disk flow and flat 
plate flow are similar in form: this is physically reasonable, because we 
would expect e€,, in general to depend on the mean distribution of turbulent 
energy (independently of any changes in sign of Reynolds stress), and this 
in turn to depend on the distribution of mean velocity profile. The velocity 
profiles are practically identical in form in the inner part of the layer, and 


so we write for the disk flow 
€_/(KU,r?) = 2%, (22) 


where the parameter p = 2g—1 = } and K isa constant to be determined, 
since this form (22) is known to give a good representation of the €,, 
distribution in the inner part of the layer in the flat plate case (5, equation 


(2)). We infer, by using (21) and (22), that in this region of the disk flow 
€_ U/(K2w2rlls) — £061, (23) 


Equation (18) can be used to show that in the inner 16 per cent of the 
boundary layer ¢% is proportional to €*7 to a high degree of approximation 
and, using (23), «,, U is seen to be proportional to Y”, which is in agree- 
ment with the result given by (20) using t = 0-54. The range of validity 
of (23) includes the region € < 0-02, where the appropriate value of the 
parameter 6 is known, and this enables the parameter A to be evaluated. 
The form y = 0-0075.X°4 can then be applied over the inner 16 per cent 
or so of the boundary layer thickness. 

We note that a theory developed on these premisses cannot be expected 
to predict accurately the distribution of temperature in the outer regions 
of the boundary layer. 

Substitution now of (20) into the temperature equation (16) leads to 


the equation amp ‘ , 
| B= al“ag) (24) 
oR Ob\" eb 
where R = r8@-95 and c = 5bwit?-5y2-95/(26—13t). Equation (24) is 
identical in form with the basic equation discussed previously by Davies 
and Bourne (5, equation (11)). We suppose that a continuous uniform 
circular source, radius ro, is situated on the surface of the disk, is concentric 
with the disk, and emits Q units of heat per unit time, per unit length of 
source. The analysis of the flat plate solution (5) can then be applied and 
a solution giving the distribution of local rate of heat transfer Q obtained 
for any prescribed radial distribution of surface temperature. In the 
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particular case of constant disk temperature, the distribution of Q(r,) is 
given by Q = K,r2, (25) 
where 

kK, (13/5)9-9'@-0(2—t)VC-OP'171 /(2—t) 'sin{a/(2—t)! 


) 
< pc, | T, —T,)b8 eS yU5, 
The total heat flux from a circular area of the disk of radius r is given by 


r 


E = | 2nr, Q(ry) dre 


sate hie 


o~ 
7 


1 > lh oF 
13 1? ze (26) 


0 
the average rate of heat transfer over a circular area of radius r by 
(arr? 10 fy. 3/5 
H = E/(ar*) = 8K, 95, 
and the Nusselt number by 


N = rH|(k(T,—T,)] = K, oR, (27) 


where 
K, 2(5/13)#2-(2—t)'@-OP 41 (2 —t) sin {a/(2—t)\bV2-9, 


and o is the Prandtl number. 

The dependence of V on o and on R is in agreement with that given 
by Cobb and Saunders from measured values; it is independent of the 
choice of 6 and ¢t. The calculated value of cA,, using 6 = 0-0075 and 
t = 0-54, is found to be 0-014, which is in good agreement with the value 
0-015 suggested by Cobb and Saunders (4) for the limiting case of al/ 
turbulent flow (Cobb and Saunders note that the accuracy to which the 
coefficient 0-015 is known is probably 5 per cent). 

At smaller angular velocities this coefficient (oA) will decrease due to 
the presence of an appreciable zone of completely laminar flow near the 
centre of the disk and also due to the presence of a significant laminar 
sublayer. An extension of the analysis is needed to include these factors, 
but the problem is complicated by the fact that the sub-layer thickness 
decreases with increasing distance from the centre of the disk (this follows 
from the definition of the sublayer thickness y, = 30p/v7,, since 7, in 
creases with increasing 7): the sublayer correction used in the flat plate 
case (5) is then difficult to apply. 

Finally, it is suggested that this method of calculating heat transfer 
in a turbulent boundary layer over a rotating disk and over a flat plate 
(previously discussed, 5) is probably capable of extension to other im 
portant problems, such as those involving a turbulent boundary layer on 
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a rotating heated sphere, although the boundary layer approximation is 
unlikely to be valid in the immediate vicinity of the equatorial plane. 
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HEAT FLOW TOWARDS A MOVING CAVITY 
By D. B. CONCER 
(Department of Physics, University of the Witwatersrand, Johannesburg)t 
[Received 20 March 1958] 
SUMMARY 

A problem arising out of the study of heat flow in deep level gold mines is that 
of the heat flow towards a stope from the surrounding rock. 

In Part I a stope is treated as an infinitely long cylindrical cavity of circular 
cross-section moving with constant velocity U, in a medium of constant uniform 
thermal conductivity A. The surface of the stope is at a uniform temperature J, 
while at points infinitely distant from the axis of the stope the temperature is equal 
to virgin rock temperature Vp, the vertical temperature gradient being neglected. 
A formula for the amount of heat entering unit length of the stope in the steady 


state per unit time is derived. 
In Part II the problem is extended to the case where the cross-section of the stope 


is elliptical. 
Part I 
1. Statement of the problem 
IN a mine, the general scheme of stoping operations is to blast away rock 
from the stope face periodically, ore being sent to the surface. For pur- 
poses of ventilation, a wall, constructed from the remainder of the rock 
not sent to the surface, is built parallel 





y to the stope face in the case of a 

parallel wall stope. In the present 

» paper the impulsive velocity of pro 

. r "gress is ‘smoothed out’ to one of con- 
tinuous motion. 

0 In this Part we will consider the 


problem of an infinitely long cylin- 
drical stope of circular cross-sectional] 
area, moving with uniform velocity U 
in a direction perpendicular to the 





generator of the stope, in an infinite 
Fic. 1, OA = radius of cavity ry. OP medium of constant uniform thermal 
aes wih Oe ae diel ney COndUCtty K, the temperature in 
of the envite. finitely distant from the axis of the 
stope being equal to J}, the virgin 
rock temperature. The surface of the stope is maintained at a constant 
temperature J). We require the steady state temperature distribution in 
+ Now at Dust and Ventilation Laboratory, Transvaal and Orange Free State Chamber 
of Mines. 


{Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 2, 1959) 
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the surrounding rock. Since the stope considered is infinitely long, the 


problem can be considered as two-dimensional, specified by the coordi- 
nates of Fig. 1. 
2. Solution 
The equation of conduction of heat to be solved is, according to Wilson 
1) 408): net7 oY rap 
(1) ev &Vv Uda 0 (1) 
dz? ' dy? ' Dox 
where D is the diffusivity of the medium, V is the temperature in the 
medium and l’ is the velocity of the cavity. (U is replaced by —U in 
the formulae quoted in the above paper.) This corresponds to a motion 
of the medium in the negative x-direction. The boundary conditions to 
be satisfied are: 
VY = J, on the surface of the cavity 2*+-y? = rj, 
Vy = JV, at infinity. 
By putting in equation (1) 
, ' Ux 
V = Vp exp/ a. 


o2 “2 2 
OF BF_ Uy _, 
cx® * Cy?) = 4 DP 


)Fex y) 


we obtain 
or, in polar coordinates, 
er 1éF 1a@F U? 
or? rar rf oP 4D 
By writing F = R(r)0(0) 
d*-) 


we obtain tn — 0 


dé? 


oR dR [Ur 
and r2 qa + r = (ais 4 n)R — 


F — 0. 


where 7 is a constant. 

Since the solution must be an even single-valued function of @, n is an 
integer, and the solution of (6) is 

Q = A, cosné. (8) 

The solution of (7) for integral n, which remains finite as r tends to 
infinity, 1s R = B,K,(2), (9) 
where K_,,(z) = K,,(z) 
is the modified Bessel function of the second kind of order n, and where 


Ur Ur 
—-z and —°=2 


2D 2D o (10) 
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The solution of (1) is then 
, , 2 cos 0 ‘. ee ee 
J Vete-* > E,K,,(z)cosné, (11) 
n x 
where the coefficients E, are to be determined from the boundary condi- 
tion at r = 1p. 
On replacing z by —iz in the expansion formulae for 


P iz cos? 


(Morse and Feshbach (2), p. 1322), we obtain 


4 


e2 008 0 = I (z)cos n@, (12) 
n is 
where f_.(z) = 2 f(z) 
is the modified Bessel function of the first kind of order n. 
Applying the boundary condition at r = r, to (11), we obtain 
(1,—V;,)exp(z, cos 4) > E, K,(z)cos nd. (13) 


Comparing coefficients of cos n@ in (12) and (13) we obtain 


ae 
E,, (I J ») - : . ( l 4) 
. , K,, (29) 
On substituting (14) into (11) we obtain, finally, 
Vy Vir | (VY, Vip) zcos 6 . To) K, (z)cos nO. (15) 


~ . K,, (29) 


n 
which is the solution of (1) subject to the boundary conditions (2). For 
large values of z, equation (15) has the asymptotic form 


v Vie ~ ¢ roost = 8 Ss T,,(20) e-*cos ne. (16) 


My Vip 22 a K,, (29) 
Now in the region $7 < 6 < 3x 
Vy 
zcos 6 Bd | 
2D 
whence (16) becomes 
4 U ix =. J (z 
ws F ~ exp| —(z— wd ti b ni 0) cos né. (17) 
cv, 2D }}\\2z K,, (29) 
n=—@ 


In the ‘wake’ left behind, 
x|+>o and y is finite. 


Uti\alt 
Therefore ih ts cae (18) 
(2D)* 
U |x Uy? 


and (19) 


~ 2D ™~ 4Diz\" 
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Substituting (18) and (19) into equation (17), we obtain 
V—V, 3 ly? 
pd f exp| —— y (20) 
where & is a function of Z,. It may be seen that (20) is an expression 
symmetrical about the z-axis, and at the centre of the ‘wake’ (y = 0), 
the temperature is VV, B 


a ae 21 
o—Ve at _ 


100° 90° goe 
110° ° r—~10° 


60° 
Soe 





Fic. 2. Variation of gro/K(Vp—V,) with 6. (There is symmetry about the line 


sis put Sued 


The amount of heat entering unit length of the cavity in unit time at 
the point (rg, @) is given by the formula 


, 


q= a. evaluated at r = rp. 
ér 


Differentiating (15) with respect to r and setting r = ry, we have 


1 lg —olexp(—29 0088) S_ Iy(zq)e08 n0— 


KU VW, —V,) . 1, (z ) ’ 
where primes denote differentiation with respect to z. 
Equation (22) may be reduced by means of the Wronskian relation 


Din(2) K(2p) — Ln (Zo) K m(Zo) = - (23) 


5092.46 Q 
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and the relationship 


ar, 
to 


m 2) 


he 1(=o) T In +1(20) 
q= oe 0) 


(24) 
exp(—Z, cos @) ‘fe ~ -. (25) 
“™ a A K,, (29) 
\ graph of the variation of qr,/AK(V;,—V,) with @ is shown in Fig. 2 for 
various values of zp. ; 
ai 








ee 


Fic. 3. Variation of the quantity Q/27K(Vp—V,) with zp. 
The total quantity of heat entering the cavity per unit length per unit 
time is given by “= 
Y qo dé (26) 
0 
, . = T(z) 
2nK(Vp—V) YS (—1)nnXe (27 
' 5 K,, (29) , 
and the variation of Q/27K(V;,,—J)) with z, is shown in Fig. 3 
Part II 
3. Introduction 


We now consider the extension of the problem of Part I to the case in 
which the cross-section of the cavity is an ellipse with its major axis 2c 
in the direction of motion; the minor axis is 2d as indicated in Fig. 4. 
4. Solution 


x 


Introducing elliptical coordinates defined by 


acoshs cost, 


y = asinhssint, 
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equation (4) becomes 
eF @F 
ot? * ds* 


where q = U*a?/16D?. 


— 2q|cosh 2s—cos 2t|F = 0, (28) 


J 


& 








0 A «x 





Fic. 4. OA ¢ = semi-major axis. OB = d = semi-minor axis. 
Motion of the cavity is in the direction Oz. 


By putting F S(s)T(t) we obtain: 


2 ’ 
: +-[2q cos 2t+-b|T = 0 (29) 
qf ** 
29 
and ; —|2qcosh 2s-+-b]S = 0, (30) 
ds* 


where / is a constant. 
Since the solution must be an even periodic function of t, the solution 
of (29) is given by 
T = A,,ce,(t,—q) (n integral), (31) 
where ce,(f,—g) is an even Mathieu function of the first kind in the 
notation of McLachlan (3). 
A solution of (30) which tends to zero as s tends to infinity is given by 


S = B, Fek,(s, —q). (32) 
The solution of (1) becomes 
V Vie 1 ¢ w cosh s cost > E, Fek,,(s, —g)ce, (t, —q), (33) 
n=—@ 
where Ua/2D = w 


and the coefficients E,, are to be determined from the boundary condition 
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at $ = 8). 


Substituting in (33) we obtain 


(Vi Vin) w cosh sq cost N E, Fek,, (89, —q)ce, (t, q). (34) 
n x 
Since [ce,,(t, —q)] 
is an orthogonal set in the range 0 < t < 27, we have 
27 
i V.—V,)f1-+3(n 
E,, o___5 [1+3(n)] | exp(w cosh 8, cos t)ce,, (t, —q) dt 
2 Fek,, (8), —q).M€. . 
0 
where 


(35) 
Ae ( ce? (t, —q) dt, 
0 
l E 0), 
and d(x) ( 


lo (4 £0). 


Following McLachlan (3) (section 2.18, equations (2) and (3)), and 
integrating term by term, we obtain the formulae 


= ¢- 
and 


( exp(w cosh 8, cos t)ce,,,(t, —q) dt = 2n(—1)" 1)" AZ”T,,(w cosh 89) 
0 — 


on 


(36) 
exp(w cosh 8, cos t)ce,,, ,,(t, —¢) dt 
0 


x 

2n(—1)" ¥ (—1) BS"4" L,+1(w cosh 89) 
r=0 

The values of Z,, are given, therefore, by the formulae 


(37) 
. m(V,—Vpp)[ 1+-8(n)] ec 
2. aw 2 1)” ¥ (—1yA2"], (w coshs 38 
2n MS, Fek,,, (89, ; q) 1) 2, 1) A? I,,(w co h o) ( ) 
and 
. (Vj —Vp)(—1)" : 
Ban +1 e : 


“wm 2n+1 


Fek, ..(8, —4) > (—1) BZ"4P Ia, 44(w cosh sy). (39) 
an+i'e r=0 > 

The temperature V at any point, (s,t), in the medium is given by (33) 

where we substitute (38) and (39) for £,,, and Ey, ,,. 

5. Flow of heat 


by 


27 


The flow of heat into the cavity per unit length per unit time is given 


gaKk {as 
J és 
0 


(40) 
evaluated at s = 8). We obtain, on differentiating (33), and placing s = s, 
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in the result, the formula 
eV 


cs 


w sinh 8, cost exp(—w cosh 8, cost) s E,, Fek,,(8, —q)ce,(t, —q)+ 
n=— © 


-exp(—w cosh s,cost)wsinhs, >} E, Fek’,(s), —q)ce,(t,—q) (41) 


where primes denote differentiation with respect to w cosh s. 
Substituting values of ce,(t,—q) given by formulae of McLachlan (3) 
(section 2.18) quoted above, and using the relation 
2costcos nt = cos(n-+-1)t-+cos(n— 1)t (42) 
and (23), we obtain 


Q = 27Kwsinhs, x 


S Ey, Feks, (8. —q)(—1)" > (—1)'A2” T',,(w cosh 8,)— 


i 
—2n7Kwsinh 8, x 


> BE, +1 Fek,,, .1(89; —q)( 1)" > (—1) BZ" 4? T5,44(w cosh 8) + 


n — & r=0 
}+- 227A w sinh 8 X 


> E,,, Fek’,,,(8, —q)(—1)" > (—1)'A2"J,,(w cosh 84) — 


n , 
—27Kw sinh 8, x 
SY Bay 4, Fekg,, .1(8, —g)(— 1)" ¥ (— 1) BRL, ..(w cosh 8). 

< r=0 (43) 
We will now verify that in the case when the ellipse tends to a circle, 
equations (33) and (43) reduce to (15) and (27) respectively. 

The limit is approached in the following way. Asa > 0, s and s, both — 00 
in such a way that 


n 


acoshs -> asinhs > r 

(44) 
a cosh 8, -> asinh 8 > ry 
and cu,,(t, —g) > cos nd, (45) 
where @ is the angle which the radius vector makes with the direction of 


motion, and AM -> (148(n))x. (46) 
Also, according to McLachlan (3) (p. 369, equations (9) and (10)) 
Fek,,, (89, —q) = Pan m—1K,,(Z9), (47) 


Fekg,, .1 (8p, oe -q) “> 8p, +1 aK, + 1(29)- (48) 
From (47) and (48) 
Feky,, (89, —q) i Pan mK, (29) ; (49) 


and Fek;,, +1(89, —q) — 89n41 71K 3n41(29); (50) 
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where in (49) and (50) primes denote differentiation with respect to the 


argument. From equations (44) to (48) we obtain 


: (Vo—Vie) Len (2) . 
£,,, > =— IK. (2 (51) 
PenT 2n (=o) 


: (V.— Ve) fen xslt : 
and | ay a n) antl 0) . (52) 
89,417 Ky, .1(Z9) 


Using equations (10), (44), (51), and (52), equation (33) becomes 


Vi.) zcos0 NS I, (0) K, (z)cos 6. (15) 


V V>+(V, 
- XH, K,(29) 


n 


Using equations (44) and (47) to (52), equation (43) becomes 


@ - 2rK(l, Vi )0| > I, (29) Lon (20) > I, 


(29) Len+1(Z0) 4 


T3,(20)Kon(Zo) Sn va(@0) Ka +10) 
a Ky, (29) . a Ky, +1(2) 


But 
> (—1)"L, (20) Ln (20) 
n x 
} 7 ( 1)"I,,(z){ L,, 1(2p) Tt I, +1(29) | 
4D (HW) (on arlZo+2 2D (—1)Ay(20)F 120) 
n x n x 
4 D> (—1)"F,(2)Ln+1(20)—4 Zz (—1)"Z,,(29)Ly +120) 
0). (54) 
Therefore 


S (yr Zal0) 7 (25) 29) + T's (29) Kul 20)} 
, 0 ( ) K ( ) n(Z9) n\2o) n(29)K,,(Zp)} 


n\~O 


Ss ( = 1)” Pn (20) Kn(20) .' 
K,, (29) 


x n x 


Ms 


( 1)" I, (29) Ly (Z0) 


ss 


T (29) Ky (29) 5 5 ( 1)" (29) 1), (29) 


> (-ye 
— K,, (Zo) 
I (29) 
— : (f, 
K,, (Zp) 


n\~O 


= - (--s)" (29) K;, (29) - Ti, (2o) Ky, (Z)}- 


(55) 
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On using equation (23), we obtain 


S 2 - 7 : . I (29) 
S (-1 I, (20) iT (2) Ka(2) + L,(z9) A, (2)} = — ~-i{P =. 
A Fe é 29) t ni o) ak ) n o) = o/s ‘ > ( ) 2 K,, (2) 


(56) 
Therefore, substituting (56) into (53) we obtain 


Q = 27K (V,,—I;) * (—1yn peo) 
Gund n\~0 


nn — <x 
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ON THE SOLUTION OF SOME AXISYMMETRIC 
BOLNDARY VALUE PROBLEMS BY MEANS OF 
INTEGRAL EQUATIONS 


I. SOME ELECTROSTATIC AND HYDRODYNAMIC 
PROBLEMS FOR A SPHERICAL CAP 


By W. D. COLLINS (King’s College, Newcastle upon Tyne) 
[Received 15 May 1958.—-Revise received 28 August 1958] 


SUMMARY 
Various boundary value problems for a spherical cap in electrostatics and in the 
theories of the potential flow of a perfect fluid and the Stokes flow of a viscous fluid 
are considered. 


1. Introduction 

In this paper we derive the solutions of certain boundary value problems 
for a spherical cap. In section 2 the electrostatic potential due to a per- 
fectly conducting spherical cap, maintained at a given potential which is 
symmetrical about the axis of the cap, is found as the real part of a certain 
integral, the expression of which is similar to that found by Green and 
Zerna (1) in the corresponding problem for a circular disk. It may be 
noted that this potential can also be found by first determining the surface 
charge density on the cap in a manner similar to that by which Copson (2) 
found the surface charge density on a circular disk. The method given in 
this paper however gives the potential directly. The potentials due to two 
given distributions on the cap are then deduced. In section 3 we show 
that the stream-function for the irrotational axisymmetric motion of a 
perfect fluid about a spherical cap, the stream-function being specified on 
the cap, can be expressed as the imaginary part of a certain integral and 
this is applied to find the stream-function for a uniform motion of fluid 
past the cap. Using an analogy derived by the author (3), we consider in 
section 4 the slow steady rotation of a cap about its axis in viscous fluid 
and find the couple required to maintain thé steady motion of the cap. 


2. The electrostatic potential problem for a spherical cap 
We require the electrostatic potential due to a thin spherical cap, which 
is maintained at a given potential symmetrical about its axis. We use 
spherical polar coordinates (r,@,¢) referred to the centre of the sphere 
r =a as origin and the axis of the cap as polar axis, the cap then being 
[Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 2, 1959] 


ON AXISYMMETRIC BOUNDARY VALUE PROBLEMS 233 


given by r=a, 0 <@< a. The electrostatic potential V(r,@) is a har- 
monic function, continuously differentiable everywhere except possibly 
at the edge of the cap, and is O(r-') at large distances from the cap. Its 
value on the cap is given as a known even function f(@) so that 


V=f(@) when r=a (O<O0< 


a). 
We now show that such a function V(r, 6) is given by 


vi.) =e { ‘ g(n)sec $n dy 


9 " "7 > 
2e'n +-a*%e-'— 2ar cos 6)* 





where g(») is a real continuous even function of 7 and where for r > a 
/(r?e*+-a?e-* — 2ar cos 6) = pet**, 
p > 0, O0O<r<a7 for 0<7 <a, 
—r7<7r<0 for —a<yn< 
and for r <a 
J (r?e! +-a?e- — 2ar cos #) = pe-*"”, 
p> 0, 0O<r<a7 for 0<7 <a, 
<0 for —-a<yn< 9. 
When r = a, we have 
, (re! +-a®e-™ — 2ar cos @) = a,(2cos 7—2co0s6) (@ > |n}), 
while 
(re + a®e-™ — 2ar cos 6) 
= +ia,/(2cos@—2cosn) (0 < |n)), 
if r > a through values greater than a, 
Fia,(2cos6—2cosy) (6 < |y}), 
if r > a through values less than a, 
(2.5) 
the upper sign holding for 0 < » <a and the lower for —a < yn < 0. 
There is thus a discontinuity in the square root at r = a (@ < |n|). We 
further note that the integral in (2.2) is to be interpreted as a Cauchy 
integral at the point r = a, 0 = 0. 

Since g() is an even function, V(r, 6) as defined by (2.2) is real and is 
O(r-') for large r. Further, the integrand of V(r,@) is a continuously 
differentiable harmonic function of r and 6 at all points where its denomi- 
nator is not zero, that is, wherever 


re +-a*e-™— 2ar cos + 0, 


this being so everywhere except on the cap r=a (0 <@ <a), hence 
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V(r,@) is harmonic and continuously differentiable everywhere except 
possibly on the cap r =a (0 <@ <a). 

We next show that V’(r, 4) is continuous for normal approach to the cap 
as r->a through values greater than a, a similar proof to the following 
holding when r -» a through values less than a. We suppose that 0 << @ < a 


and note the relation 


r2et aze~* 2arcos 6 2 

$a*r*(cos 7 —cos 0)? + (7 a)*|(r a)*--4ar(1—cos 6 cosy) | 

ta*r*(cos ” cos @)?--(r—a)?. 
Thus for r > a we have 

r-e Late Par cos 6 » Pa? cos ”) cost. (2.6) 
and hence obtain 
q(n)see 4n q(n) see hy 
(r*e') +-a*e-") — Zar cos 0) ay 2c0s 7 —2 cos 6 ; 


By the theorem of dominated convergence we have that, as r—->a 
through values less than a, 
x x 
q(y)see by dy F ; gi n )sec } n dy (2.7) 
J y(r’e'+-a*e-™— 2arcosé) J a,(2cos 7—2cos@) 
x —@ 
provided g(y) is bounded, this holding at all points on the cap except 
# = 0. This failure at the point r = a,@ = 0, does not however invalidate 
the solution, since it is shown later that V(r, 0) tends to a finite limit at 
this point. The square root on the right of (2.7) is to be interpreted as 
the limit of that on the left in accordance with (2.3) and (2.5). 
We further require that V(r, @) be continuously differentiable on either 
side of the cap.t The zeros of the function 
reetS+arze-%—Jarcosé (¢ n +1€), 
regarded as a function of ¢, lie in the upper half {-plane when r > a and in 
the lower half plane when r < a. We can therefore displace the path of 
integration in (2.2) from —a to «a away from the real axis into the lower 
and upper half ¢-planes for r > a and r <a respectively, provided we 
assume g(C) to be an analytic function in a simply-connected domain which 
contains both the old and the new paths of integration. In the new regions 
of definition of V(r, #) the integrands are continuously differentiable func- 
tions of r and 6, other than on the edge of the cap r = a, 0 = a. Therefore 
the derivatives of V(r, @) with respect to r and @ exist and are continuous 
in the neighbourhood of the cap, and can be calculated by differentiation 


+ I am grateful to a referee for certain suggestions in regard to this and the preceding 
parts of the proof. 
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under the integral sign. Thus V(r,@) is harmonic in this neighbourhood 
and tends to finite limits as the point (r,@) approaches a point on the cap, 
including the point (a,0), in any manner through values of r greater or 
less than a, these limits being equal by the argument of the preceding 
paragraph. Since V(r,@) admits this harmonic extension through the cap, 
it must in its original region of definition be continuously differentiable 
on either side of the cap, except at its edge. 
Using (2.3) and (2.5), we have from (2.7) on r a(O <@< a), 

4 

r —-g(»)see 4y dy 


V(a,0 7 : - {(@). 
J 4 (2cos 7—2 cos @) Je) 


0 


Making the substitutions ¢ = tan}, s = tan $6, and writing 
q() G(tan $n), f(?) F(tan 38), 


we reduce this equation to 


wk oe (0 <8 < tan $a). (2.9) 


[ Gi(t) dt F(s) 
(et) (+8?) 


This equation is now of the form considered by Green and Zerna (1) and, 
if f(#) is continuous for 0 < 6 < «, its solution is 


2d sF (8) ds 


Git) 
‘ ” dt j ((f@—s —s?)(14 -s)}! 


(0 <t < tan $a). (2.10) 


When G(t) has been determined from (2.10), the value of V(r, 6) follows 
from (2.2). 

When a = 7, the cap becomes a sphere of radius a. The expression 
(2.2) still holds for either the region r > a or the region r < a. It can be 
shown from (2.10) that, if /(@) is continuously differentiable for 0 < 6 < z, 
g(n) is O(cos* }n) for 7—« < » < 7, where e is a small positive quantity. 
It then follows that V(r,@) is continuously differentiable at all points for 
which r > a (or r <a) including the point r = a, @ = z. 

The total charge & on the cap is the limit, as r —- 0, of (4) 





ira | g(n)sec 4n dy 


(r2e*7 +-a2e-* — 2ar cos 6)! 


so that 2=a [ 9) dy. (2.11) 
0 


As illustrations of these results we consider the following problems. 
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(a) Spherical cap at constant potential. For a spherical cap at a constant 
potential U’ we have f(@) = U and from (2.10) we obtain 


Ot) = ———.., 
\ a(1+t?) 
2U 
or g(y) = —cos? $y. 


7 


On substituting in (2.2) and integrating, we obtain 


U 
V(r, 0) = —(ry+ay’), (2.12) 
nr 
siny siny 2sina 
where - = = Sa caapeemaes 
a r rit+fe 
and r? = r?+a*—2arcos(a—8), rz = r?+-a?—2ar cos(a+8). 


The angle } y(r,@) is such that 0 < y < 4m at all points (r,6,d4) other 
BIC Y Y ) I 


than those which lie in that segment of the sphere r = a bounded by the 
surface of the cap and the plane, rcos @ == acos a, containing its edge, for 
which points 7 < y < 7. Further, if S is the sphere 


r? cos a—ar cos 0 0, 


through the edge of the cap and the centre of the sphere r = a, the angle 


va y'(r,@) is such that 4a < y’ < a for those points (7,6,¢4) exterior to 
) y 2 ) } 


9 


the sphere r = a but interior to S for « < 47 and exterior to S for a > 4}. 
For all other points 0 < y’ ; har. 
From (2.11) the total charge = on the cap is given as 
. Ua ; 
a (x-+-SIN a); 


so the capacity C of the cap is given by 


0 >» a(a+sin «) 


U 7 





(b) Spherical cap in a uniform field. If an earthed spherical cap is placed 
in a uniform field £ parallel to the axis of the cap, we write 


V = —Ercosé-4 V, 
so that on the cap the perturbation part V, of V satisfies the condition 


V,= Eacosé@ (0< 0 <a). 





9 
8° 


Hence F(s) = Fa(—*}, 
and from (2.10) we find that 
2 Fa(1—3t2) 


OY 4 Peppa, oe all 
) a(1+t?)? 


b] 
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2Ea 
so that y(n) = ae cos }y cos $y. 


7 


On substituting in (2.2) and integrating, we obtain 





T ” ee r3 


where z = rcos@ and the angles y and y’ are interpreted as in (2.12). 
These results are in agreement with those obtained by alternative 
methods by Ferrers (5), Gallop (6), Basset (7), and Szegé (8). 


, ae azy’  a*(r? cos a—az)tan y’ 
V(r, @) B:+=|271 (a COS x —z)tan y +——~ _ > y ’ 


3. The stream-function for the motion of a perfect fluid past a 
spherical cap 
In an irrotational axisymmetrical motion of a perfect fluid the Stokes 
stream-function ¥ satisfies the equation 
ce sindée/ il @ 
BOY es Boe hee een ee ee @ (3.1 
or? r® @0\sin 6 06 
together with the condition that is a constant over any rigid boundary 
present in the flow. 
We now suppose that the stream-function ‘Y(r,@) for such a flow past 
a spherical cap is given by 
V(r, 0) = bo(r, 0)+(r, 4), 


where u(r, @) is the stream-function for the flow in an infinite unbounded 
fluid in the absence of the cap and ¢(r, @) the perturbation stream -function, 
which must be added to y,(r,@) to make the boundary of the cap, r = a 
(0 < @ <a), a stream-sheet. We require to determine y(r,@) at a general 
point of the fluid. On the boundary of the cap ‘Y’(r, @) is constant and thus 
d(r,@) is a given even function of 6, say f(#), on this boundary. We have 
W(r,0) = f(0) on r=a (0 <6 <a), (3.2) 
where /(@) = constant —y,(a,@). Since the constant is at our disposal, we 
choose it so that f(@) = 0 at the point r = a, 0 = 0. 
Since the velocity components q,, gg, in the r- and 6-directions, given by 
1 oF 1 oF 


q7,.=- --=o d= : SS 
dr r2sin@ 6° 16 rsin@ ér’ 





are continuous at all points except possibly those on the edge of the cap, 
u(r, @) must be continuously differentiable at these points. Further, since 
the introduction of the spherical cap into the flow can cause no additional 
total flow over the sphere at infinity, the velocity components derived 
from the perturbation stream-function are O(r-*) at a large distance r from 
the origin and hence y(r, @) is O(r-') at this distance. 
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We now show that a stream-function u(r, @), satisfying these conditions. 
is given by 
x x 


* a(n) (re? cos 6— ae 11) dn r * q(y)(ae "cos 8—rei?”) dn 
Pids(r, @) : _ ‘- # 


ane in Yar cos b)! a (r2e” j ane a Par cos @) 
. 
(3.3) 


where g(m) is a real continuous odd function of 7 and the square root 


(re 


+ a%e-™— 2ar cos 6)! is interpreted as in (2.3), (2.4), and (2.5). 

Since g(7) is an odd function of », the stream-function y(r,#) given by 
(3.3) is a real function. The expression A* can be evaluated at all points 
where the necessary differentiations can be carried out under the integral 
signs, that is, at all points where the integrands so obtained are continuous 
functions of (7,4, +). The only points where this condition is not satisfied 
are points such that 


r2e+™ + a2e*' — Jar cos 6 0, 
that is, points for which r = a, y -@. All such points lie on the cap 
} a(oO<@ x). Now the integrand of the first integral is readily seen 


to be a solution of (3.1) and the application of Butler’s inversion theorem 
for a sphere (9) to this integrand shows that the integrand of the second 
integral is also a solution of (3.1). Hence ¢(r, @) is a continuously differen 
tiable function of r and @ and satisfies (3.1) everywhere except possibly on 
the cap r a (0 G x). 

As for the electrostatic potential in section 2, we can show that ¢(r, @) 
is continuously differentiable on either side of the cap except at its edge 
r a,@ x, and that the limits of the integrals in (3.3) as (r, #) approaches 
a point on the cap are equal to the values of the integrals at that point. 

We next note that, if the integrands in (3.3) be expanded as power 
series in a/r, we obtain 


asin’? ,. fa\® 
yi(r, A) g(y)sin n(1+-cos n) dy + terms in | } , (3.4) 
r r 


showing that y/(r, #) has the required order for large r. Using (2.5), we have 
on the cap 
))sin 4 
ws(a, @) 2(1-+-cos @) gi7)sin 37 on f(®) (r=a,0 <6 <a). 
J (2 cos 7—2 cos 6)! 
0 
Making the substitutions ¢ = tan }n, s = tan $0, and writing 


g() G(tan 4»), {(@) F(tan 36), 
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as before. we obtain 


7 tGi(t) dt _ wr ; ‘ 
a. - 1 F(s),(1+s8?) (0 <8 < tan 4a), 3.5 
(+f) (ee) ¢ (ayia) ? An 2) (3.5) 


0 


which is again an integral equation of the type considered by Green and 
Zerna (1). If f(@) is continuous, its solution is 


t 
tG(t d [i -+ 8") F(s) d: 
) 1S (QUteVOe 6 <;. tan }a). (3.6) 


lt? 2ndt : , (B—s?) 
) 





When G(f) has been determined from this equation, ¢(r,@) follows from 
(3.3). We note that, if F(s) is given as a constant, G(t) is O(t-*) for small ¢. 
Hence it is necessary to choose the disposable constant in F(s) = f(@) to 
ensure that F(s) = 0 at the point r=a, 6= 0. We then have F(s) 
is O(s*) for small s and it follows from (3.6) that G(t) is bounded at the 
origin. 

The cap becomes a sphere when a = 7 and the expression (3.3) then 
applies to motions external or internal to the sphere. In this case it 
follows from (3.6) that, if f(@) is O(cos? 46) for 7—« < @ < 7, where e is 
a small positive quantity, g(y) is O(cos}y) for 7—d < » < 7, where 6 is 
a small positive quantity depending on e. It can then be shown that 
U(r, #) is continuously differentiable at all points for which r > a (orr < a) 
including the point r = a, @ = 7. Now the stream-function (7, @) for the 
motion in the infinite unbounded fluid is constant between any two sources 
on the axis of the motion, the value of the constant changing by an amount 
2m at such a source of strength m. For motions external to the sphere 
any such singularities lie in the region r > a and thus f(@) is zero at r = a, 
# = m, since it is zeroatr = a,@ = 0. It then follows that f(@) is O(cos* 40) 
near @ = 7. For motions internal to the sphere the total change in y(r, 4) 
along that part of the axis inside the sphere is 2x (total strength of the 
sources on the axis inside the sphere) and this must be zero since fluid is 
conserved inside the sphere. It then follows that f(@) is zero at r = a, 
6 = 7, and is thus O(cos* 36) near 6 = z. 

We next give a corresponding result to the above for a circular disk of 
radius ¢ with its centre on the axis of symmetry. If (w, ¢,z) are cylindrical 
polar coordinates with the centre of the disk as origin and the axis of 
symmetry as z-axis, we require a stream-function ¢%(w,z) such that on 
the disk p=f(w) (7: =0,0<w<e), 
where f(a) is a continuous even function of w with f(0) = 0. Further, 
U(aw,z) is a solution of (3.1), is continuously differentiable at all points 
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except those on the rim of the disk z = 0, w = c, and is O(r-!) at large 
distances r from the origin. It can be shown by a similar proof to that 
for the cap that (aw, z) is given by 

. 

2iub(w, 2) | 


. 
c 


g(t)(z+-at) dt 


ar? t-(2+-at)2te’ 
where g(t) is a real continuous odd function of t, such that 


t 
wm. 2 i ee 
. a dt a @ (t? —a?) 

0 
The square root {w?+-(z-+-it)?}! is interpreted similarly to the square root 

occurring in (2.2). 

We now give an illustration of the use of (3.3) and (3.6): 

A spherical cap in a uniform stream. If a spherical cap is placed in a 
uniform stream of velocity U in the negative z-direction, the axis of the 
cap being taken as the z-axis, the perturbation stream-function u(r, 4) 


satisfies F 


wb = f(0) kU a*sin?é (r a,0 <@ x). 
at : 2U a*s? 
Thus fk (Ss) 20? 
(1+8*)* 
and from (3.6) we find that 
2U at 
G(t) one 
a(1-+4-¢*) 
Var. — 
or g(n) = — -sin 7. (3.7) 
On using (3.3) we obtain 
U ca oe a 
u(r, @) = s yr* sin?@— y’ — sin*6+-(r cos 6+-a)(a cos x—r cos 6)tan y+ 
aT r 


a , 
+ —(r-+acos @)(r cos a—a cos @)tan y | (3.8) 
> 


where r,, 72, y, and y’ have the same meanings as in (2.12) 


4. The slow steady rotation of a spherical cap in a viscous fluid 

If a solid of revolution rotates slowly about its axis with angular 
velocity © in a viscous fluid, the only non-zero velocity component is the 
component v in the ¢-direction. The author showed in a recent paper (3) 
that solutions for v are given by 


ub 


: 4.1 
rsin@ ( 
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where y is the stream-function for the same solid moving with constant 
velocity U along its axis in perfect fluid, provided U is replaced by —2Q 
in #. This result can now be applied to give the solution for a rotating 
spherical cap. In fact (3.8), besides giving the perturbation stream- 
function for a uniform motion past the cap, also gives the stream-function 
for a translation of the cap with uniform velocity U in the positive z- 
direction. Replacing U by —2Q and using (4.1), we can thus derive v 
from (3.8). 

In order to find the couple Q required to maintain the rotation of the 
cap we use (3.4) to find that, at a large distance r from the origin, v has 
the form ‘ 

v= = | g()sin 7(1+-cos n) dy, 
0 
where U’ has been replaced by —2Q in the function g(y). If u is the vis- 
cosity of the fluid, the couple Q required to maintain the motion is given 
by (10) 


Q = —2r7p lim sintort = (" dé = 87pa | g()sin n(1 +-cos ») dy. 
—> 2 r\r 
: 0 
Hence from (3.7) we have 
a 
Q = 16Qya { sin?y(1+-cos n) dy = 4Qua5(2a—sin 2a + $sin*a). 
0 
If « tends to zero and a tends to infinity in such a way that aa tends to a 
finite limit c, the couple Q becomes that required to maintain the rotation 
of a circular disk of radius c, namely, 
32Qc3 
9 = ane 
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ON THE EQUATION OF A SYNCHRONOUS MOTOR 


By W. A. COPPEL (Birkbeck College, London) 


[Received 10 July 1958] 


SUMMARY 
The strongly non-linear differential equation #+ a#+-sinx x7, Which arises in 
the study of synchronous motors, is discussed for small values of «a by the method 
of small parameters and by a generalization of the method of slowly varying ampli- 
tude and phase. 


1. Introduction 
THE oscillations of a self-starting synchronous motor are governed by the 
non-linear differential equation (1) 


¢+ae%+sinxe = B, (1) 
where « and £ are positive constants. The same equation also represents 
the motion of a pendulum which is acted on by a constant torque and 
opposed by viscous friction. The first mathematical study of the equation 
(1) is due to Tricomi (2), some of whose results are reported in Andronov 
and Chaikin (3) and other texts on non-linear mechanics. Among more 
recent investigations those of Amerio (4) and Bohm (5) may be especially 
noted. Amerio (6) has also generalized many properties of the equation 
(1) to a wide class of equations of the form 

“4 F(x, x), 
where F is periodic in 2. 

The question of greatest interest is the asymptotic behaviour of the 
solutions when ¢—» +00 and, in particular, the manner in which this 
behaviour depends on the initial conditions and on the values of the 
constants a, 8. This question has been completely answered in a qualita- 
tive sense by the work of Tricomi and Amerio. It has been shown that 
there are two possible types of solution—the ‘convergent’ solutions which 
tend to a finite limit as t-> +0, and the ‘divergent’ solutions which 
behave asymptotically as the sum of a linear function and a periodic 
function. Divergent solutions exist if and only if the first-order equation 

y dy/dx+ay+sinax B, (2) 
which is obtained from (1) by putting y = #, has a positive solution y(z) 
of period 27. If B > 1 all solutions are divergent. If 8 < | and if the 
ratio 7 = B/« is less than a critical value »,, depending on £, then all 
solutions are convergent; if 8 < 1 and » > », then both convergent and 
divergent solutions are present. The initial conditions which characterize 
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the different types of solution can be determined at once from a knowledge 
of the two solutions of (2) which pass through the singular point 
(7—sin~'f, 0). 

On the quantitative side all that is known concerns the critical value 
n(8). It has been calculated to three decimal places for B = 0(-1)1 by 
Urabe (7) and Giger (8), by numerical integration of the first-order equa- 
tion (2). The difficulty in obtaining further information has been attributed 
to the strongly non-linear character of the equation (1). In fact, since x 
is unrestricted, it is not permissible to replace sinx by x or x—}25, as 
one often does on other occasions. 

In the technical applications, however, « is usually small and this 
suggests that we can treat (1) not as a quasi-linear equation, but as a quasi- 
conservative one. This is what is done here. It is shown that the rigorous 
method of small parameters of Poincaré (9) for the determination of 
periodic solutions can be applied to (2), and provides new information 
about the divergent solutions.+ It is shown also that the heuristic method 
of slowly varying amplitude and phase, which has been developed by 
Kryloff and Bogoliuboff (11) for certain quasi-linear equations, can be 
extended to the equation (1), and leads to equations soluble in closed 
form. The results obtained in this way are in agreement with the known 
qualitative results. 


2. The method of small parameters 
2.1. Consider first of all the general equation 
y dy/dx+-g(x) = af(x,y), (3) 
where g(x) is a continuous function of period w such that 


[ g(x) dx = 0, 

a 
and where f(x,y) is a continuous function with period w in x having a con- 
tinuous partial derivative with respect to y. We wish to know if, for small 
values of «, this equation has a solution y(x) which is positive and of 
period w. 

For «a = 0, (3) reduces to the elementary equation 
Y dY /dx+-g(x) = 0. 
This has the positive periodic solutions 
Y = Y(z,C) = 24|C—G(z)}}, 


where G(x) = (ea) dé (4) 
0 


+ The method of small parameters has already been used by Urabe (10), but in quite 
a different manner and not for the determination of the periodic solutions. 
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and C is any constant greater than max G(x) (0 < x <w). Since Y takes 
the value (2C)! at x = 0, we have C = 4$Y*(0). 

Suppose now that a + 0, and let ¥(x) be a solution of (3) which is defined 
and positive for 0 <2 <w. For y(x) to form part of a solution of period 
w it is necessary that y(w) y(0), and since the coefficients of (3) have 
period w this condition is also sufficient. But the differential equation (3) 
is equivalent to the integral equation 


by?(x)— by?(0) + G(x) x | LE. med) dé, 
0 


and therefore, since G(w) 0, y(w) y(0) if and only if 

w 

| flax, y(x)| da 0), (5) 

0 
Thus (5) is a necessary and sufficient condition for y(x) to have period w. 

We now ask for what values of C does there exist a solution y = y(2, «) 

of (3), defined for all sufficiently small «a and having period w in x, such 
that y(x,a) converges to Y(x,C) uniformly in 0 <2 <w asa—>0? By 
making «->0 in the relation (5) we see that © must be a root of the 
equation a 


$(C) = [ flx, ¥(x, C)] dx = 0. (6) 
0 


Conversely, suppose C = C* is a root of this equation and write 

Y *(x) Y (zx, C*). 
For all sufficiently small a and for all y, sufficiently near Y *(0) the solution 
y = y(x, a, ¥,) of (3) which takes the value y, at « = 0 will be defined and 
positive for 0 <2 <w. We will obtain periodic solutions converging to 
Y*(x) if we can choose y, = y (x) so that the integral 


w@ 


ub(a, Yo) S [x y(x, x, Yq) | dx 
0 
vanishes and so that y)(a) > Y*(0) as a+0. By the implicit function 


theorem this will be possible in one and only one way if 
-_ 
py (0, ¥*(0)) A 0. 
Since W(0, Yo) = (492) 


and @Y /éy, = y,/Y we have 


Ww 


ie » fifa, ¥*(x)} 
b, (0, ¥*(0)) = Y*e pte ’ 
by ( (0)) (9) Yaz) dx 
i 
Hence, if d(C) = 0 and if - 
"a [x, Y(x, C)| 
€ = bh a a? L 
si } Y(x,C) mes 


0 
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then for all sufficiently small values of « there is one and only one solution 
of (3) of period w in the neighbourhood of Y(z, C), and this solution con- 
verges uniformly to Y(x,C) as «a > 0. 

We consider next the stability of such a periodic solution. If we put 
z = y*, (3) is replaced by the equation 


dz/dx = F(x, 2). 
where F(x,z) = —2g(x)+-2af(x, z*) 


has period w in x. It is known—see, for example (12)—that a periodic 
solution 2(2) of an equation of this form is asymptotically stable or un- 
stable for x -» +00 according as the equation of first variation 

dl/dx = Fix, 2(x)\f 
has negative or positive characteristic exponent, that is, eecording as 

a 

| Fjx,z(x)| dx $0. 

0 
Returning to the original variable y it follows that for all sufficiently small 
a the periodic solution found above will be asymptotically stable or un- 
stable for 2 —- --co according as a¢'(C) is negative or positive. 


2.2. From now on it will be supposed that f, is everywhere negative and 
that for each fixed value of x there is a value of y for which f(x, y) is negative. 
These are essentially the assumptions made by Amerio (6). The first 
condition implies that f is a strictly decreasing function of y. Hence, by 
the theorem of monotonic convergence, ¢(C) is a continuous, strictly 
decreasing function of C for which —a < ¢(%) < 0. 

We shall show that if ¢(C) < 0 then for all sufficiently small a all 
solutions of (3) of period w lie in the region y < Y(x,C). In fact we can 
choose C’ < € so that also d(C’) < 0. Let y(x) = y(x, a) be the solution 
of (3) which takes the value (C+-C’)! at x = 0. For all sufficiently small 
: this solution will be defined for 0 < 2 < wand will satisfy the inequalities 


¥ (a, C’) < y(x) < ¥(z,C). 
Therefore, since f is a decreasing function of y, 


w 


| fix, y(x)| dx < d(C’) < 0. 


0 


If 7(x) is a solution of (3) of period w, it must always be greater or always 
less than y(x), since the continuity of f, ensures that only one solution 
passes through any point in the half-plane y > 0. But it cannot always be 
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greater than y(zx), since this would imply 
w w 
| fl[x,j(x)| dx < | flx, y(x)| dx < 0, 
0 0 
which contradicts (5). Therefore 
(x) < y(x) < Y(z,C). 
Similarly it may be shown that if 4(C) > 0 then for all sufficiently small « 
all solutions of (3) of period w lie in the region y > Y(a, C). 

Let M=maxG(x) (0<24<w). (7) 
If d(M) > 0 then the equation 4(C) = 0 has a unique root C = C* greater 
than M, and ¢’(C*) < 0. Since for 5 small and positive 

¢(C*—3) > 0 > g(C*+8) 
any solution of period w must lie between Y(x, C*—8) and Y(x, C*+8). 
Moreover, by the results already established, there is exactly one solution 
of period w in this region for all sufficiently small «. 

On the other hand, if ¢(M) < 0 the equation ¢(C) = 0 has no root 
greater than M. In this case (3) has no positive solution of period w, for 
all sufficiently small «. In fact any positive solution of period w must lie 
in the region y < Y(x,C), for every C > M, and hence in the region 
y < Y(x,M). But since Y(x,M) certainly vanishes at some point this is 
impossible. 


The results which have been established can be summed up as follows: 


If d(M) > 0 then for all sufficiently small «, (3) has a unique positive solu- 
tion of period w. It is asymptotically stable or unstable for x > +0 accord- 
ing as «20. As «1-0 it converges uniformly to Y(x,C*), where C* is the 
unique root greater than M of the equation d(C) = 0. 

If 6(M) < 0 then for all sufficiently small «, (3) has no positive solution 
of period w. 


Let us now apply these results to the equation of a synchronous motor. 

The differential equation (2) will be written in the form 
ydy/dx+-sinx = a(n—y), (2)’ 
where a and 7» = f/« are positive constants. Thus in the present case 
g(x) = sin, f(z.y) = n-y, 

and the various assumptions which we have made are all satisfied. 
Evidently G(x) = 1—cosx = 2sin? iz, 
and hence : 


d(C) 277 - -(2C)! | | 1—(2 ‘C)sin® ha}? dx 
0 


2nn—Sk-1 E(k), 


ON THE EQUATION OF A SYNCHRONOUS MOTOR 247 


where k = (2/C)* and E(k) is the complete elliptic integral of the second 


kind. Moreover M = max Gz) = 2, 


and ¢(M) = 22yn—8. 
Substituting these values in the theorem above we arrive at once at this 
conclusion: 

If » < 4/m then for all sufficiently small a the equation (2)’ has no 
positive solution of period 27, but if 7 > 4/7 then for all sufficiently small 
« it has one and only one such solution. This solution is asymptotically 
stable for x + +00 and converges uniformly as a > 0 to 

(2/ky)(1—k§ sin? $2)!, 
where k, is the unique root between 0 and | of the equation E(k) = }ank. 

The more complicated equation obeyed by a salient pole synchronous 
motor can be dealt with in exactly the same way. 

The solution 2 = x(t) of the original equation (1) corresponding to any 
solution y = y(x) of the phase plane image (2) is obtained by inverting 
the relationship 


J dé 


-t — : 
, y(€) 


In particular, to a positive periodic solution of (2) corresponds a solution 


of (1) which is the sum of a periodic function and a linear function. There 
is no difficulty in obtaining an approximate representation of this solution, 
using the limiting form of the periodic solution y(x) as «-—> 0. However, 
this will not be carried out since the result will later be obtained by a 
different meth 


2.3. To obtain a closer approximation to the periodic solution of (3) 
we follow the procedure of Linstedt (9). Suppose, for simplicity, that 
f(x,y) is an analytic function of y. Then the periodic solution y = y(z, a) 
is analytic in « and we can write 

Y = Yot ay, +a7y2t+..., 
where the coefficients y; all have period w. It follows that 
Y= Yor ayy totyet., 
af (x,y) = af (x, Yo) +f, (2, Yo) +-.- 
Substituting these expansions in (3) and equating coefficients of like 
powers of « we get 
YoYo+9(x) = 0, 
YUTMYo = f(&, Yo); 
Yost WnVtYe2¥o = Lyle Yol¥r» 


7 . . . . 
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The first equation gives us 
Yo = 24 cq -G(x)]! 
for some constant cy. From the second equation we get 


x 


Yo" = 17 ( SLE, wolE)] dé 


0 


for some constant c,. Since y, must have period w we arrive, as before, 
at the equation 


@ 

{ fl. yolx)] dx = 0 
0 

to determine cy. The third equation gives us 


YoYot dy? = C2+ | FALE vol) lal) dé 
0 


for some constant c,. Since y, must have period w we obtain the equation 
w 
| Fil x. Yo(x) lys(2) dx = 0 
0 
to determine c,. And so on. 
For the equation of a synchronous motor w = 27 and f, = —1. Hence 
the last equation can be written 
rf 
| y,(2) dx = 0. (*) 
wade 
The equations for y, and y, become 
Yo = (2/ky)(1—K sin? $2)}, 
“3 
YoY = Cy + na—(2/ky) (1—k@ sin? $€)! dé. 
0 
Since y, is an even function the condition (*) reduces to c, = 0. Thus, 
neglecting powers of « above the first, the periodic solution of (2)’ is given 
by zr 


y = (2/k,)A(x)- x{ inky x - | A(é) dé} /A(z), 
0 


where k, has its former meaning and A(x) = (1—kgsin? }z)!. 


3. The method of averaging 

3.1. The method described in the previous section has the disadvantage 
that it provides information only about the periodic solutions. To deduce 
the behaviour of all the solutions the qualitative results of Tricomi and 
Amerio must be used. We are going to apply now a different method 
which is not subject to this restriction, but which has no pretence to 
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rigour. It is a generalization of the method which Kryloff and Bogoliuboff 
(11)—-following van der Pol—have applied to the quasi-linear equation 

+a = f(x, #). (8) 
For purposes of comparison we will first briefly describe their method. 
The equation ee 
obtained from (8) by putting « = 0, has the solutions x = acos(t+94,). 
where a and @, are arbitrary constants. Kryloff and Bogoliuboff therefore 
introduce new variables a and & defined by the equations 
x= acosy, 
z= —asiny. 
It is then easily shown that the second-order equation (8) is replaced by 
the pair of first-order equations 
ad = —esing f(acosy%, —asin’s) = «F(a, h) 
w = 1—ea cosy f(a cosy, —asinys) = 14+-€G(a, sb) J 
The functions F(a, %) and G(a, ys) are periodic in % with period 27. They 
can therefore be broken up (uniquely) into two parts 


(9) 


ce 


F(a, $) = F(a) += Fila. ), 


b 


G(a, 4) = Gla) +5 G,(a,), 


where 


2a 

, l al 

Fla) = 5 | Flaw) dy, 
0 


2a 


l 

G,(a) = — [ Gia, wb) db 

Zr | 
0 

are the mean values of F and G respectively, and where F,, G, are periodic 

in & with period 27 and have zero mean value. Suppose now that « is 

small. Then the equations (9) show that formally 


a=Oc), w= 14+O0(e). 
Therefore ee 
os 


Thus the equations (9) can be rewritten in the form 


a = €F,(a)+-« 4 Fy(a, ys) + Oe?) 
( 


j= Lh eGila) +4 Gla, W) +02) 
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The approximation procedure consists first of all in neglecting the terms 
of order e? in these equations. Secondly, the average trends a and J of a 
and % are obtained by neglecting the small periodic terms involving F, 
and G,. Thus @ and ¢ are defined as the solution of the system of equations 
a = «F,(a), 
2 oe 
yb =e. i+ e(7, (a). 
These values are then substituted for a and y% in the right-hand side of 
(9)’. The resulting equations can be integrated at once and give 
a = 4+ eF,(4,d), 
ub yt e(7,(d, pb). 
These are the complete first-order approximations. However, for many 
purposes the motion is already adequately described by the average trends 


dand wd. The effect of the second terms is simply to add small oscillations 
about this mean motion. 


3.2. We return now to the equation of a synchronous motor and write 
it in a form analogous to (8), namely, 


x+sinx uf (%), (10) 
where f(z) = n-2z. (11) 


To put the equation in a more suitable form for approximating we apply 
the method of variation of constants to the well-known solution, in terms 
of elliptic functions, of the equation for an ordinary pendulum: 
¢+sinz = 0, 
Let W = W(t) = sin? r+ }2?; (12) 
W is constant for the pendulum equation, 2W being its total energy apart 
from an arbitrary additive constant. We take different forms for the 
solution of (10) according as W 2 1, corresponding to the two types of 
pendulum motion—oscillation and revolution. 
if W > 1 and & > 0 we write the solution in the form 


sinjr=snu, cos}r=cnu, }¢ = k-'dnu, (13) 


where k and wu are functions of the independent variable t, and k is the 
modulus of the elliptic functions. For this change of variables to be a 


legitimate one the transformation must be invertible. But from the 
definition of W we have 


W = sn?'u+k"dn?u = k-?, (14) 
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and consequently the inverse transformation is given explicitly by the 


equations k = (sin® 424+ }2%)-}, 


iz 
“= | (1—k* sin2£)-! dé. 
0 
From (14) we get by differentiation 
d 


a) = }2(2sin }rcos ha +-2). 
€ 


By (10) the right-hand side is equal to }a#f(#). Hence 


ok) = ak" dnu f(2k-'dnu). 


Also, differentiating the first equation (13), we get 
d . 
7 u) = sécos $x = k-enudnu. 
¢ 


But by the formulae for the total differentiation of elliptic functions (13) 


d du sg d 
— (s = —-—4Sdu. - 
(snu) = enudn ee 3 Sdu 7 


k?)}, 
dt (©) 


where, in Neville’s notation, 


u 
2 
Sdu = | ae M du 


dn? u 
0 


Combining this with the previous equation we get 
du 


d 
oe «x k-14 1 Bde. 2). 16 
di k-1+438Sdu a? (16) 


Thus the second-order equation (10) in the variable x is equivalent to the 
two first-order equations (15) and (16) in the variables k and w. 
The functions F(k,u) = kdnuf(2k-"dnu), 
G(k,u) = Sdu—(u/K)Sd K 
are periodic in wu with period 4K, where K = K(k) is the complete elliptic 
integral of the first kind. The same also holds for their product. More- 


over, since /'G is an odd function of u its mean value is zero. 
Suppose now that a is small. From the equations (15) and (16) we have 


d p29) _ du 4 
na ) = O(a), ye k-*+- O(a). 


Proceeding as in the method of slowly varying amplitude and phase, and 
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using the fact that Sd K = 2dK/d(k*) (13), it follows that 
> oF (k)+a . F(k, u) + O(a), 
tf 
u dK 
K dt 
2K 


[ k-" dn u f(2k-! dn u) du 


~2K 


l S G,(k, u) + O(a2). 


l 
Here F(k) = — 

, 4K 
is the mean value of F, while F, and G, are periodic in « with period 4K 
and mean value zero. To obtain the average trends & and @ of k and u 
we neglect the terms of order e? and the small periodic terms involving F, 
and G,. Thus k and @ satisfy the equations 


i. . . in 
= () = aF,(k), (17) 


du = K-14 wd dk (18) 
dt K dt 

Having obtained k and @ the complete first-order approximations for 
k and u can immediately be found. However, since the complete expres- 
sions are complicated and add nothing essential, in what follows we shall 
restrict the discussion to the average trends k and a. Moreover, for con- 
venience of writing we will drop the bars and write simply k and u. 

The solution of the first-order linear differential equation (18) is 


t ’ 
oe Ir 
= K(t a. 19 
, | k(r)K (7) ais 


7). For the case of a synchronous motor, that 


Consider next equation (1 
is for f(z) given by (11), 


] 
F(k&) = — 
1(4) ik 


> 

| k-' dn u(n—2k-! dn u) du, 
—2K 
or, since dnu = dn(—u) = dn(2K—u), 

K 
F,(k) = z | k-' dn u(n—2k-! dn u) du. 
0 

The integral here is easily evaluated—see, for example, (14)—and we 
find that 


ae) = 2a(tarnk—E)/k*K, (20) 


where E = E(k) is the complete elliptic integral of the second kind. This 
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can be written in the form 
3k-3Kd(k?) 
tan—kk 
But using the formulae for the derivatives of complete elliptic integrals 
we have d(k-1E)/d(k®) = —}k-3K. 
Hence the solution of (20) is given by 
k7E = }rn+ce,e™, (21) 


= —adt. 


where c, is a constant. 

Now E£ decreases when k increases and takes the value 1 at k = 1. 
Hence if » < 4/7, the function E—}7nk is always positive. In this case 
it follows from (20) that W = k-? — as ¢t increases, and from (21) 
that after a finite time we must have W = 

On the other hand if » > 4/7 the ee E = }rnk has a unique root 
k, between 0 and 1, and d(k-*)/dt 2 0 according as k 2 ky. Hence k con- 
verges monotonically to ky, whether initially it is greater or less than this 
value. By (21) the time taken to reach the limit k, is infinite. Moreover, 
writing (21) in the form 

darn —4kg 2K o(k?—k2) + O(k—k,)? = }en+c,e-, 
we see that for t > +00 
k—k, = 8, e~at 
elles 0 
Similarly for t > +0 
t 


-2 
+e ) HaK@ de Kw f Fed 


= (t— to) /kyt-(At+ B)e-™, 
where ¢, is a constant of integration and A and B are constants whose 
values can easily be determined. 
3.3. These solutions are for the case W>1,z2>0. If W>1 and 
% < 0 the form of the solution must be modified. We then take 
sinjz = —snu, cosjx=cnu, $4 = —k-"dnu. 


Equations (14) and (16) remain unchanged, but (15) is replaced by 


5 (k-*) = —ak-!dnu f(—2k-'dnw). 


Consequently the equation (20) for & becomes 
© (kt) = —2a(benk+ BYR. 


The solution of this equation is obtained by changing the sign of 7 in (21), 
kB = —}nn+e,e™. 
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It follows that W — k-2 decreases as ¢t increases and that after a finite 
time we have W l. 


3.4. Suppose finally that W < 1. In this case we write the solution in 
the form . , 
f sinjr=ksnu, cos}xu=dnu, 4$¢ = kenu, (22) 
where again k and u are functions of t, and k is the modulus of the elliptic 
functions. For definiteness it will be supposed that initially —7 < x < z. 


From the definition of W we have 


W = k*sn?'u+k?en?u = k?. (23) 
Proceeding as for the case W > 1 we find that 
d (k?) = akenu f(2kenu), 
dt 
du fc sn u a 
> 1+4{Sd «— ——_—___—__ } - — (K°). 
di " k?cnudn | di\ 
By splitting up the right-hand sides in the same way as before we get 
” pee , ” e 
atk) uF (k) : 7 Fy(k, u)+O(«?), 
du u dk d 
- 1+ - - ++ G, k, a}. 2 . 
dt K dt “di a “)s a ) 
2K 
] 
*here F(k) ome ( kenu f(2kenu) du 
where 1( iK | enuf ( ) 
-~2K 


and F,, G, are periodic in uw with period 4K and mean value zero. For 
the equation of a synchronous motor, that is for f(z) given by (11), 


2K K 
Coe 2 
F,(k) ate kenu(n—2kenu)du= —— ( k? en? u du 
: 4K 3 l K J 
2K 0 


2(E—k’?K)/K, 


where k’* = 1—k*. Thus the average trends k and @ satisfy the equations 


(k) —2a(E—k’*K)/K, (24) 
du u dK 
a ee ; 25 
dt K dt _ 


The solution of the first-order linear differential equation (25) is 


t 
dr 


p= Ke [ ay (26) 
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Also, from the formulae for the derivatives of complete elliptic integrals, 
d(E—k'®?K)/d(k?) = 4K. 
Hence the solution of the separated variables equation (24) is given by 
E—k'*K = c,e-%, (27) 
where c, is a constant. 
Now the function E—k’?K increases from 0 to 1 as k increases from 
0 to 1. Therefore, by (24), W = k* decreases as ¢t increases, and by (27) 
k? + 0 ast—> +oo. Since |sin}x! < k < 1 it follows that x remains inside 
the interval (—7z,7) and tends to zero as t-» +o. (In general 2 will 
converge to the nearest even multiple of 7.) When ¢ is large we have 


from (27) 
oe Po 


7 


ke = 


Similarly from the —— 


-. om 
-1 


we get u=t—t, ri. — (t—ty + a-!)e- 


where /, is a constant of loci 

It only remains to piece together the various solutions which have been 
obtained. Suppose that a solution for which W > 1 initially arrives at 
W=1 when t=t,. There are three possibilities consistent with the 
differential equations (20) and (24) which W satisfies in the regions 
W 21: 

(i) W passes immediately into the domain W < 1, 

(ii) W remains equal to 1 for all t > t,, 

(iii) W remains equal to 1 for an interval t, <t < ¢, and then passes 
into the domain W < 1. 

The approximate equations naturally do not enable one to decide 

between these possibilities. However, by appealing to the exact equation 
(1) we can show that W cannot be constantly equal to | throughout an 
interval of time, and on this basis we can exclude cases (ii) and (iii). Hence 
our conclusions may be summarized thus: 
If » < 4/7 all solutions converge to a finite limit as t > +-00, but if n > 4/7 
only those solutions converge for which initially  < 0 or sin? 42+-}a@* < 1. 
The remaining solutions behave asymptotically as though x were the angle 
made with the downward vertical by a simple pendulum of length g which 
performs complete revolutions in a time 2k, K(ky), where ky is the unique 
root between 0 and | of the equation E(k) = 4ank.t 


+ Since k, decreases as » increases the period of revolution is a monotonic decreasing 
function of 7. 
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‘These results agree with the qualitative results stated in the introduc- 
tion. There are two types of solution, the convergent solutions which tend 
to a finite limit as t-» +-o, and the divergent solutions which behave 
asymptotically as the sum of a linear function and a periodic function. 
Moreover there is a critical value », of the ratio 7 = B/« such that all 
solutions are convergent if and only if 7 < »,. The limit of the convergent 
solutions is given as zero (mod 27), whereas the exact value is easily shown 
to be sin-'8 (4). Thus the approximations give the limit of the exact value 
for 8 +0. The initial conditions which characterize the divergent solu- 
tions when 7 n, may also be regarded as the limit for 8 + 0 of the 
known qualitative results. Finally, as proved in section 2, the value 4/7 
for the critical value », and the value 2k,K, for the time period associated 
with the divergent solutions are the limits for 8 + 0 (B/a = constant) of 
the exact values of these quantities. 
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SUMMARY 


By using a Karman—Pohlhausen method the distribution of local rate of heat 
transfer is evaluated for the case of air flow in an axisymmetrical laminar boundary 
layer ona heated circular cylinder, the temperature of the cylinder being independent 
of downstream distance. This calculation serves to link the numerical values 
obtained by Seban, Bond, and Kelly for small downstream distances to those ob- 
tained by Bourne and Davies for large downstream distances. 


1. IN a previous paper by Bourne and Davies (1) the influence of curvature 
on the calculated values of local rate of heat transfer from a heated surface 
was demonstrated by considering the problem of forced convection from 
a heated circular cylinder into the surrounding, incompressible, axisym- 
metric laminar boundary layer produced by a uniform stream. An exact 
solution of the temperature equation had already been given by Seban 
and Bond (2), with a correction by Kelly (3), in the form of a series, valid 
for small values of va/Ua® (< about 0-04), where U denotes uniform main- 
stream velocity, 2 downstream distance, v kinematic viscosity, and a radius 
of the cylinder. Bourne and Davies obtained an asymptotic series solution, 
which is valid at large values of vx/Ua? (> about 100) and is based on the 
type of analysis given by Glauert and Lighthill (4). 

In order to bridge the gap, for the intermediate values of vx/Ua?, 
between the results given by these two series solutions, Bourne and Davies 
applied in a modified form an approximate method of calculation which 
they had previously introduced for the case of heat transfer from a flat 
plate. This method was based on a power law approximation to the 
Pohlhausen velocity profile and, although capable of application to the 
general case of arbitrary distribution of cylinder temperature, suffers from 
the disadvantage of requiring a separate power law fit for each point on 
the bridge. The maximum deviation of the power law from the appropriate 
Pohlhausen profile (over 95 per cent of the boundary layer thickness) 


(Quart. Journ. Mech. and Applied Math., Vol. XII, Pt. 2, 1959) 
5002.46 Ss 
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increases from about 4 per cent at the smaller intermediate values of 
vz/Ua® to about 10 per cent at the larger values. In the latter region the 
power law method is unlikely to lead to really accurate results.t 

A check on the results obtained by this method is now achieved in this 
paper by an alternative calculation. A Karman-—Pohlhausen approxima 
tion is employed, following the method described by Glauert and Lighthill 
in calculating the distribution of skin-friction on the cylinder, and the 
distribution of heat transfer obtained by numerical integration of a first- 
order linear differential equation. This method is likely to be more accurate 
at the larger intermediate values of vax/Ua?. 


2. The temperature integral relation, expressing conservation of heat 
between two planes perpendicular to the axis of the cylinder, is 
4 - 
pC, - u(T'—T,)(a+y) dy = (S) a, (1) 
n y=0 


0 
where 7' is the temperature of the air stream at normal distance y from 
the surface of the cylinder, 7, the temperature of the main stream, u the 
downstream velocity, 5,, the thickness of the temperature boundary layer 
at downstream distance 2, p the fluid density, c,, the specific heat of the 
fluid at constant pressure, k the thermal conductivity of the fluid. 

The discussion given by Glauert and Lighthill (4) shows that their 
Pohlhausen velocity form represents flow conditions near the cylinder 
surface fairly accurately. This may be written in the form 

u ] ‘ 
—=-In(l+y/a) for y <6 = a(e*—1) 
l : : (2) 
and u/U=1 for y>5 
where U is the mainstream velocity, 5 is the thickness of the momentum 
boundary layer at downstream distance z, and a is a function of va/Ua? 
given numerically by Glauert and Lighthill. Their argument is imme- 
diately applicable in the context of heat flow to show that the temperature 
distribution is well represented by the forms 
T7,—T 1 
—l___ — —|n(l+y/a) for y <8, = a(e®—1) 
7,—T,~ B si 

T,—T 


where 7; is the temperature of the cylinder surface and 8 is a function of 
va/Ua® to be determined, corresponding to « in equation (2). Substituting 


+ An accuracy of about 5 per cent was suggested by Bourne and Davies (1). 
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(2) and ” into (1) leads to the relation 


4 | lt, sin(l Ly a) [inc +-y/a)\(a+y) dy + 
8a 
+] ! —in(1-+yla) (a y) dy\ = 
re) B | 


4 
te 
for 5,, > 4 (corresponding to values of o, the Prandtl number, < 1), where 
« is the thermal diffusivity of the fluid. A similar relation can be obtained 
for 5,, < 6 (corresponding to o > 1). Writing y = a(e?—1) in equation 
(4) and integrating yields 


ak fe22{ 48(—4) —Ma?—a-+-4)]+4(1-+8)} + 


Differentiating and using the relation 


x-*{(2a2 <= 4, 
dx Ua? 


given by Glauert and Lighthill, we find finally that 


iB 1(g2a e28) } 228 __ at al -1)}] 
x 


2(20-!— 1)e?*+- (14+ 8— —ao~!){a-? —¢e2%(q-2—2a-1)}. (7) 
This equation is now in standard form and is to be solved subject to the 
end condition 8 = 0 at a = 0, and at large downstream distances the 
values of 8 must run smoothly into the values given by the asymptotic 
solution. This enables 8 to be calculated as a function of «, which is 
already given numerically as a function of vz/Ua*® by Glauert and Lighthill. 
The local rate of heat transfer Q, per unit length of the cylinder, is then 
given by 2rk(T,—7,)/B. 


3. For small values of va/Ua?, we note that « is given in terms of a 
series in the variable vx/Ua? by Glauert and Lighthill (4), and 8 in terms 
of another series in (va/Ua?) by Seban and Bond (3). In this region the 
values of « and £ can thus be calculated for selected values of va/Ua? and 
the values of 8 corresponding to, say, « = 0-05, 0-10, and 0-15 can be 
obtained by interpolation. With these initial values and the condition 
a == 0,8 = 0, equation (7) can be integrated numerically using a standard 
step-by-step forward integration procedure. In this paper the Milne- 
Simpson method has been used, the steps in « being 0-05 initially and 0-1 
for values of « greater than 0-3, and the integration taken as far as a = 3-0 
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(corresponding to log,,(va/Ua*) = 2, in the range of the asymptotic series 
solution). The results shown in Table | were obtained by interpolation, 
using the relation between «a and va/l’a® given by Glauert and Lighthill. 


TABLE | 


Calculated values of local heat transfer rates in air, per unit length 
of a circular cylinder, obtained by two methods 


ote aa 
eis k(T, —T,) 














Values obtained by numerical integration of Values obtained by the power 
( va | the (Karman—Pohlhausen) differential equa- law method of Bourne and 
? lia tion (7) Davies 
1-0 i 0-97 0-96 
0-5 0-82 0-79 
0 0-68 0-64 
Od 0-56 0-52 
1-0 | 0-46 0-42 
| 








The numerical results for Q/k(7,—7,), shown in Table 1 over the inter- 
mediate range of values of va/l’a? (corresponding to the bridge), are found 
to run smoothly into the results given by the asymptotic series solution 
for log,,va/Ua? > 2, and are seen to be in fairly close agreement with the 
results obtained by the power law method employed previously by Bourne 
and Davies (1). As noted above, the error involved in representing the 
Pohlhausen profile by a power law increases, as vax/U’a? increases, and the 
power law method, when applied at the downstream end of the bridge, is 
therefore not likely to lead to very good accuracy. On the other hand, 
the results given in the second column of Table 1, by numerical integra- 
tion of (7), are of the same order of accuracy as those given for skin 
friction by Glauert and Lighthill, and we suggest they should replace the 
previously calculated values given in column 3 of the table. 
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NOTE ON ALLEN AND SOUTHWELL’S PAPER 
“RELAXATION METHODS APPLIED TO DETERMINE 
THE MOTION, IN TWO DIMENSIONS, OF A VISCOUS 

FLUID PAST A FIXED CYLINDER’ 


By MITUTOSI KAWAGUTI 


(Institute of Science and Technology, University of Tokyo, Japan) 


SUMMARY 


In the present note, the results of Allen and Southwell’s paper, which dealt with 
the motion of a viscous fluid past a fixed circular cylinder by applying relaxation 
techniques, have been compared with the results hitherto known, both theoretical 
and experimental, and it has thus been shown that there are some discrepancies 
between Allen and Southwell’s results and the other authors’. It seems that the 
origin of these discrepancies may be found in the unsatisfactory treatment of the 
point at infinity or in the use of coarse meshes in Allen and Southwell’s study. 


1. Recentriy Allen and Southwell (1) have dealt with the two-dimen- 
sional steady flow of a viscous fluid past a circular cylinder, by the use 
of the relaxation method, They have thus obtained numerical solutions 
for various values of Reynolds number, R = 0, 1, 10, 100, 1000, assuming 
the flow to be laminar. We can infer the characteristic features of their 
solutions from many figures in their paper, although the details of their 
computation were omitted, e.g. the treatment of the point at infinity, the 
mesh length, the accuracy of their solutions, ete. 

There are both theoretical and experimental papers dealing with the 
same problem, especially for the lower Reynolds number régime where 
the flow is stationary. Thus, Thom (R = 10, 20) (2) and the present 
author (R = 40) (3) obtained the numerical solutions for the same problem 
by the use of Liebmann’s method. As the starting form of such a numerical 
solution, the present author adopted the results derived from his study on 
the problem by means of the Galerkin method (4). The author also studied 
the limiting form of the laminar solution when the Reynolds number 
becomes infinitely large (5). For the very low Reynolds number régime, 
Tomotika and Aoi (6, 7, 8) treated the same problem using Oseen’s 
approximation. On the other hand, Thom (2), Kovasznay (9), Homann 
(10), and Taneda (11) studied the flow experimentally. 

The characteristic features of Allen and Southwell’s solutions seem to 
be inconsistent with the results hitherto known: 
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(i) In Fig. 1 are shown the values of the characteristic length a/r of 
the standing vortices behind the cylinder (cf. Fig. 2) which have been 
obtained from Allen and Southwell’s paper, in comparison with the other 
results. It will easily be seen that Allen and Southwell’s solutions show 


a quite different feature from the others’. According to the author’s 





a 
=m " 
y a 
A 
2 Ux 
a 
_—--@--—-—- oe 
a ee G 
& = | 
- 
a | 
Oe . ee . L pepeet 
| 10 100 R 1000 
Fic. 1. The characteristic length of the standing vortices behind the circular 


cylinder. 
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opinion, the ratio a/r should become infinitely large as the Reynolds 


number increases indefinitely (3, 4). 

(ii) In Fig. 7, Allen and Southwell gave the pressure distributions 
over the surface of the cylinder. These pressure distributions have their 
maxima at points other than the forward stagnation point where the 


pressure distributions obtained by the others show their maxima. 
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In view of these discrepancies, the author cannot but be suspicious of 
the accuracy of their results. The origin of the discrepancies between 
Allen and Southwell’s and the other authors’ results may be found in the 
unsatisfactory treatment of the point at infinity or in the use of coarse 
meshes in Allen and Southwell’s study. 


REFERENCES 


1. D. N. pe G. ALLEN and R. V. SourHweEL Lt, Quart. J. Mech. App. Math. 8 (1955) 
129. 
2. A. Tuom, Proc. Roy. Soc. A, 141 (1933) 651. 


3. M. Kawacutt, J. Phys. Soc. Japan 8 (1953) 747. 

4. Report Inst. Sci. Technol., Univ. of Tokyo 2 (1948) 33 (in Japanese). 

5. J. Phys. Soc. Japan 8 (1953) 403. 

6. S. Tomorika and T. Aor, Quart. J. Mech. App. Math. 3 (1950) 140. 

¥ — Memoirs Coll. Sci., Univ. of Kyoto A, 26 (1950) 9. 

8. - Quart. J. Mech. App. Math. 4 (1951) 401. 

9. L. S. G. KovAsznay, Proc. Roy. Soc. A, 198 (1949) 174. 

10. F. Homann, Forschung auf dem Gebiete des Ingenieurwesens 7 (1936) 1. 
11. S. Tanepa, J. Phys. Soc. Japan 11 (1956) 302. 








CORRIGENDA 


L. 8. D. Morey, ‘An improvement on Donnell’s approximation for 


Vol. XII, Pt. 1, 1959. 


thin-walled circular cylinders’, 
References 3 and 4 refer to the J. Appl. Mech. 


Liv, ‘On the separation of gas mixtures by suction of the thermal 


i 
Vol. XIT, Pt. 1, 1959. 


diffusion boundary layer’. 


Page 11, footnote. The word ‘temperature.’ should be added to complet 


the sentence. 
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